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FOREWORD 


This  technical  report  was  prepared  for  the  Air  Force  Rocket  Propulsion  . 
laboratory,  Edwards  Air  Force  Base,  California,  by  Rocketdyne,  a 
Division  of  North  American  Aviation,  Inc.  The  report  covers  Phase  II 
of  Contract  AF04(6ll)-9714  in  accordance  with  Part  I,  paragraph  B.2.b.(3) 
under  Project  No.  3058,  Task  No.  305802,  and  Program  Structure  No.  750G. 
The  Air  Force  Project  Monitor  was  Lt.  L,  A.  Barlock,  office  symbol  RPREE. 
The  research  reported  was  conducted  during  the  period  1  June  1964  through 
31  May  1965  by  the  Research  Department,  R.  J.  Thompson,  Jr.,  Vice  Pres¬ 
ident  and  Director.  The  Rocketdyne  Program  Manager  was  R.  B.  Lawhead, 

\ 

withv J.  M.  Zimmerman  and  T.  A.  Coultas  as  responsible  supervisors. 

J.  D.  'Seader  was  technically  responsible  for  the  overall  effort. 

H.  A.  Friedman  was  the  Principal  Investigator,  an'd  with  B.  L.  McFarland 
developed  the  mathematical  procedures  utilized  in  the  resulting  digital 
computer  program.  S.  Persselin  and  the  Principal  Investigator  were 
responsible  for  coding  the  program. 

This  report  has  been  given  the  Rocketdyne  identification  number  R-6050-2, 

This  report  contains  no  classified  information  extracted  from  other 
classified  documents. 

Publication  of  this  report  does  not  constitute  Air  Force  approval  of  the 
report,  findings,  or  conclusions.  It  is  published  only  for  the  exchange 
and  stimulation  of  ideas. 


ABSTRACT 


An  alternating  direction  implicit  finite  difference  procedure  is 
developed  for  solving  a  general  class  of  two-dimensional  transient 
non-linear  ablation  an*,  heat  conduction  problems  for  rocket  engine 
thrust  chamber  walls.  Engine  firing  may  be  steady  or  can  consist  of 
multiple  start  or  pulsing  type  modes.  The  analysis  is  performed  in 
cylindrical  coordinates  (axial  and  radial)  assuming  an  axisymmetric 
multi-material  wall  having  arbitrarily  curved  boundaries  and  interfaces. 
A  maximum  (due  only  to  computer  storage  limitations)  of  five  different 
wall  materials  may  be  treated  in  a  given  problem  including  one  charring 
ablative  material.  Temperature  dependent  properties  may  be  specified 
for  each  material.  Chemical  reactions  within  the  ablative  material 
are  handled  in  depth;  also  surface  erosion  due  to  chemical  reactions 
or  melting  at  the  hot  gas  boundary  is  treated  and  the  resulting  surface 
recession  is  predicted.  The  numerical  procedures  have  been  programmed 
in  Fortran  IV  for  automatic  computation  on  the  IBM  7094  digital  com¬ 
puter.  Comparison  of  the  results  of  sample  computations  with  actual 
engine  test  filing  data  is  included  in  this  report. 
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INTRODUCTION  AND  SUMMARY 


In  order  to  effectively  design  rocket  engine  thrust  chambers  which 
are  partly  or  completely  cooled  by  heat  sink  or  ablative  techniques 
or  evaluate  engine  firing  data  from  such  chambers#  it  is  necessary 
to  adequately  predict  the  effect  of  combustion  gas  products  upon  the 
wall  materials  for  the  particular  geometrical  configurations  of  interest. 
In  1963,  in  conjunction  with  Contract  No.  AF04(61l)-8190  (Ref.  1  ), 
Rocketdyne  formulated  a  transient  one-dimensional  mathematical  model 
of  the  mechanisms  affecting  charring-ablative  wall  materials,  i.e., 
heat  transfer,  erosion,  and  corrosion.  However,  a  detailed  solution 
of  the  equations  was  not  attempted  at  that  time. 

In  April  1964,  Rocketdyne  initiated  a  12-month  two-phase  program 
entitled  "Effect  of  Rocket  Engine  Combustion  on  Chamber  Materials" 
under  Contract  AF04(61l)-9714»  In  Phase  I  of  this  program,  which  is 
reported  in  detail  in  Ref.  2  ,  two  goals  were  achieved.  First  of  all, 
the  equations  of  the  1S6J  model  were  verified.  Secondly,  a  computer 
code  based  on  an  augmented  non-linear  transient  ono-dimensional  model 
was  written  to  relate  the  combustion  environment  to  its  effect  on 

1 

ablative  thrust  chamber  materials.  ' 

In  Phase  II,  as  described  in  detail  in  this  report,  finite  difference 
procedure  and  accompanying  computer  program  was  developed  for  solving 


a  general  two-dimensional  extension  of  the  Phase  I  problem.  The  analysis 
was  performed  in  cylindrical  coordinates  (axial  and  radial)  assuming  an 
uxisymmetric  multi-material  solid  with  temperature  dependent  properties 
and  arbitrarily  curved  boundaries  and  interfaces.  A  maximum  of  five 
wall  materials  can  be  handled  in  the  computer  program  including  no  more 
than  one  charring  ablative  material.  Surface  erosion  of  the  materials 
due  to  chemical  reactions  or  melting  at  the  hot  gas  boundary  is  also 
treated  in  tho  program  and  the  resulting  surface  recession  is  predicted. 
All  heat  fluxes  at  the  boundaries  and  material  interfaces  are  expressed 
in  terms  of  normally  directed  temperature  gradients. 

There  are  already  in  existence  several  one-dimensional  and  two-dimensional 
programs  for  analyzing  ablation  and  heat  conduction  problems,  but  none 
it  is  felt  provides  adequate  methods  for  handling  these  normal  gradient 
conditions  at  curved  boundaries  and  interfaces.  Furthermore,  in  most 

i 

previously  existing  programs,  an  explicit  forward  difference  analog  is 
used  to  simulate  the  energy  equation  in  each  time  step.  In  order  to 
avoid  numerical  instability  with  this  method,  it  is  often  necessary  to 
employ  prohibitively  small  time  steps.  In  the  present  program,  on  the 
other  hand,  a  generalization  of  the  unconditionally  stable  implicit 
Peaceman-Rachford  alternating  direction  method  (see  Ref.  3  for  basic 
technique)  has  been  employed  to  discretize  the  energy  equation  in  con¬ 
junction  with  special  second-order  accurate  techniques  for  simulating 
normal  gradient  conditions  at  curved  surfaces.  Coupled  to  the  solution 


of  the  energy  equation  for  the  temperature  distribution  in  each  time 
step  is  a  finite  difference  solution  to  the  first  order  t vo-dimensional 
mass  continuity  equation  for  the  mass  flow  rate  of  the  generated  pyrolysic 
gases  in  the  internally  charring  ablative  material*  The  gas  mass  flow  is 
assumed  to  be  zero  at  each  interface  separating  the  region  from  other 
solid  materials* 

In  formulating  the  model  and  developing  the  numerical  solution,  the 

i 

following  assumptions  were  made  concerning  the  charring  ablative  materials 

1.  There  is  assumed  to  be  a  single  continuous  charring  material 
which  is  described  by  a  continuous  rather  than  an  interface  model.  Thus, 
the  moving  pyrolysis  zone,  which  separates  the  virgin  material  from  the 
char  residue,  is  distinguishable  only  by  the  temperature  range  within 
vhich  the  pyrolysis  reaction  takes  place. 

2.  Chemical  reactions  occurring  in  the  charring  region  are  treated 
in  depth  but  are  simulated  thermodynamically  in  the  heat  storage  and 
gas  enthalpy  terms  of  the  energy  equation.  They  include  a  maximum  of 
three  gas  generation  reactions  plus  cracking  of  the  generated  gases.  In 
addition  to  pyrolysis,  a  c-har-reinforcement  reaction  and  a.  further  decom¬ 
position  of  the  solid  product  of  this  reaction  can  also  be  handled. 

3.  The  porous  char  is  cooled  convectively  by  the  pyrolysis  gases, 
the  gas  and  char  temperatures  at  any  point  being  assumed  identical. 
Conduction  of  heat  within  the  gas  is  taken  to  be  negligible  in  comparison 
with  that  in  the  char. 
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4.  All  momentum  effects  and  thus  consideration  of  a  momentum 
equation  are  neglected  due  to  the  two  following  asaumptionsi 

a)  negligible  gas  density  relative  to  that  of  the  solid  material 

b)  orientation  of  the  gas  mass  flow  vector  in  the  char  with 
the  temperature  gradient  rather  than  the  pressure  gradient' 
vector  (under  certain  conditions,  further  restrictions  are  - 
placed  on  the  direction  of  flow). 


The  heating  and  cooling  mechanisms  incorporated  in  the  model  and  handled 
by  the  program  include  engine  firing,  environmental  heating,  radiation  to 
the  environment,  insulation,  chemical  reaction  and  melting  of  the  hot 
gas  wall  surface,  blocking  of  the  hot  gas  by  the  generated  and  tran¬ 
spiring  gases  at  the  boundaries,  chemical  reaction  and  convective  cooling 
(discussed  above)  in  the  charring  material  and  conduction  within  caeh 
wall  material  and  across  each  material  interface.  Intermittent  engine 
firing  with  soakback  is  also  treated  including  both  steady  firing  and 
high  frequency  pulsing.  The  latter  is  simulated  in  representative  time 
periods  by  reduced  values  of  the  convective  heat  transfer  coefficient. 

The  principal  parameters  computed  and  printed  out  at  preassigned  time 
levels  are  the  temperature  distribution  throughout  the  wall  materials, 
the  gas  mass  flow  rate  at  the  hot  surface  of  the  chairing  ablator,  and 
the  new  surface  positions  of  the  receding  hot  gas  boundary.  Also  computed 
and  optionally  printed  out  are  effective  heat  transfer  coefficients  at  the 


trail  boundaries  and  (both  axial  and  radial)  mass  flow  rates  within  the 
region  of  charring.  Finally,  a  matrix  la  supplied  at  the  termination 
of  the  calculation  of  the  maximum  temperature  achieved  at  each  compu¬ 
tational  point  within  the  charring  material.  This  is  particularly  use¬ 
ful  in  case  of  multiple  start  or  pulsing  operations.  Input  requirements 
include  the  temperature  dependence  of  the  physical  and  chemical  properties 
of  the  wall  materials,  the  axial  variation  of  the  various  surface  heat 
flux'  and  hot  gas  surface  erosion  parameters,  and  explicit  piecewise 
quadratic  functions  describing  the  nozzle  wall  configuration.  Also 
required  are  the  initial  thermal  conditions  as  well  as  the  composition 

of  the  charring  material. 

, 

/ 

z' 

Program  options  includes  (l)  up  to  5  wall  materials,  only  one  of  which 
may  char  internally,  (2)  erosion  of  the  hot  gas  surface  of  any  of  the 
materials  exposed  to  fixing,  (3)  omission  of  bilker  charring  or  eroding 
materials  which  results  in  a  straight  conduction  analyais. 

In  this  report  a  brief  description  is  given  of  previous  work  in  the 
area  of  two-dimensional* transient  ablation  calculations  and  a  comparison 
is  drawn  with  the  present  program,  a  discussion  and  derivation  are 
presented  of  the  ablation’ and  heat  conduction  equations,  the  numerical 
procedures  developed  to  solve  the  equations  are  described,  and  the 
essentials  of  the  computer  program  are  discussed.  The  results  of  several 
program  computations  are  presented  and  compared  with  available  experimental 
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test  data,  and  typical  exampj.es /Of  computer  runs  are  described  in  detail- 

r 

Finally,  recommendations  .are  given  for  further  development  .and  imple¬ 
mentation  of  the  two-dimensional  program. 


PREVIOUS  INVESTIGATIONS 


The  widespread  utilization  of  ablative  materials  for  nose  cones  of  re-r 
entry  vehicles  and  for  rocket  engine  thrust  chambers  has  stimulated  the 
development  of  analytical  methods  for  predicting  the  transient  behavior 
of  ouch  materials.  As  described  in  Ref-  2  ,  a  number  of  analytical 
methods  have  been  developed  f  or  handling  one-dimensional  models  in 
either  rectangular  or  radial  coordinates. 


On  the  other  hand,  relatively  little  effort  has  gone  into  the  develop¬ 
ment  of  multi-dimensional  models  which  are  needed  for  a  more  accurate 
description  of  heat  transfer  phenomena  in  the  complex  geometry  of  actual 
multi-material  thrust  chamber  walls.  In  what  follows  we  give  a  brief 
summary  of  previously  developed  multi-dimensional  models,  as  discussed 
in  available  literature,  including  the  use  of  thermal  analyzer  computer 
programs,  models  developed  primarily  for  nose-cone  applications,  and  a 
recant  two-dimensional  nozzle  model.  A  comparison  is  then  drawn  of  the 
capabilities  of  the  computer  code  developed  in  the  present  Phase  U 
program  with  the  others  described.  -  - 

THERMS!  ANALYZES  COMPUTES  PHO  SEAMS 

Eocketdyne  has  made  extensive  use  of  modified  versions  of  two  thermal 
analyzer  digital  computer  programs,  namely  TIGER  (Hef .  4  )  and  TAP 
(Ref.  5  ).  Both  of  these  programs  permit  the  solution  of  multi- 


dimensional  heat  conduction  problems  for  multi-material  walls  and 
*  ' 

structures.  Surface  boundary  conditions  may  consist  of  radiation  and/or 

convective  me ’as  of  heat  transfer.  An  explicit  finite  difference  solution 

* 

technique,  subject  to  the  usual  stability  criteria,  is  utilized  in  both 

.  / 

/ 

programs. 

/ 

The  TAP  program  is  the  more  versatile  of 'the  two  since  it  can  handle 
temperature-dependent  physical  properties,  a  wider  .variety  of  surface 
boundary  conditions,  and  can  more  easily  approximate  cylindrical  co¬ 
ordinate  systems.  However,  like  the  TI&SR  ■program,  it  is.  not  suf- 

t  * 

ficiently  versatile  to  solve  the  problems  ,of  interest  'to  the  present 
program,  namely  those  involving  solid  materials  which  are  subject  to 
internal  charring  and/or  surface  ablation,  without  considerably  over¬ 
simplifying  the  actual  physical  and  chemical  processes.  Furthermore, 
the  user  of  the  TAP  program  must,  after  selecting  the  positions  of  the 
nodes,  submit  as  input  to  the  program  the  coordinates  of  each  node  as 
well  as  all  the  corre  oonding  capacitances  and  conductances.  For  complex 
geometries,  the  necessary  pre-calculations  as  well  as  the  writing  of 
data  sheets  can  be  quite  time-consuming.  As  will  b6  discussed  in  detail 
later,  the  computer  program  presented  in  this  report  automatically  sets 
up  the  node  points  so  as  to  minimize  the  required  pre-calculations. 

MULTI-DIMENSIONAL  COMPUTATION  METHODS  FOR  NOSE  CONE  TYPE  PROBLEMS 
Hurwicz,  et.al.  (Ref.  6  and  7  )  have  described  multi-dimension  computation 
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schemes  for  electronic  analog  computation  which  have  been  developed 
to  solve  particular  simultaneous  ablation  and  heat  conduction  problems 
for  re-entry  structures.  In  order  to  handle  surface  recession,  a 
shrinking  coordinate  system  is  utilized  in  conjunction  with  the  usual 
finite  difference  equations.  Typical  calculations  for  several  structures 
indicate  the  significant  differences  between  more  exact  multi-dimension 
calculations  and  one-dimensional  approximations. 

TWO-DIMENSIONAL  NOZZLE  MODEL? 

McCuen,  et.al.  (Ref.  8  )  have  developed  a  two-dimensional  single¬ 
material  computer  program  primarily  for  application  to  non-melting, 
non-charring  nozzle  throat  inserts.  The  surface  exposed  to  the  hot 
combustion  gases  can,  however,  recede  due  to  chemical  reactions  con¬ 
trolled  either  by  mass  transfer  or  reaction  kinetics.  A  moving  co¬ 
ordinate  system  is  used  in  conjunction  with  an  explicit  finite  differ¬ 
ence  techniques..  Only  one  wall  material  may  be  specified,  but  a  variety 
of  internal  wall  boundary  conditions  are  possible.  Temperature-dependent 

material  properties  can  be  utilized  including  directional  thermal 

> 

conductivities  to  simulate  anisotropic  materials. 

\ 

COMPARISON  WITH  PRESENT  PROGRAM 

The  numerical  procedures  described  below  and  the  resulting  computer 
code  developed  in  Phase  II  of  the  present  program  are  more  comprehensive 


o 
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in  many  important  respects  than  those  previously  available  programs 
discussed  above.  Charring,  for  example,  is  handled  in  the  present 
program  as  a  t*n>-dimensional  continuous  process  and  not  a3  a  dis¬ 
continuous  interface  reaction;  on  the  other  hand,  except  for  the 
Hurwicz  program,  which  is  not  applicable  to  rocket  nozzles,  charring 
is  omitted  in  the  other  programs.  In  this  connection  it  should  be 
mentioned  that  a  version  of  TAP  has  been  adapted  at  Rocketdyne  to 
simulate  charring  by  effectively  lumping  all  gas  generation  and  con¬ 
vection  effects  in  the  storage  term  of  the  energy  equation.  Of  the 
programs  discussed  above,  only  the  McCuen  code  includes  the  effect  of 
boundary  ctirvature  in  the  normal  temperature  gradient  conditions  at 
the  hot  gas  surface,  by  all  means  the  most  important  region,  through 
use  of  a  radial  curvilinear  transformation.  The  present  code  can 
account  for  curvature  at  all  boundaries.  The  present  program  employs 
an  implicit  '  -nitc  difference  procedure  to  solve  the  energy  equation 
thus  avo  'o’  ‘miting  the  size  of  the  time  step  used  in  order  to 
insure  the  stability  of  the  numerical  computation.  This  limitation 
exists  in  the  other  programs  available  except  for  that  of  Hurwicz 
which  requires  no  discretization  of  the  time  derivative  at  all,  being 
an  analog  procedure.  Usually  an  implicit  method  requires  a  time  con¬ 
suming  iterative  solution  in  each  time  step,  but,  as  discussed  below, 
the  alternating  direction  method  (ADM)  used  in  the  present  program 
generates  a  tridiagonal  system  of  equations  in  each  time  step  which  can 
be  solved  directly  by  well  known  recursion  formulas.  Although  the 


explicit  method  is  between  two  and  three  times  faster  than  the  ADM  in 
each  time  step,  the  ADM  can  often  use  a  t  .me  step  as  much  as  100  times 
larger  than  that  needed  for  the  explicit  method  and  is  therefore  general¬ 
ly  much  more  efficient.  One  excellent  fc  'ture  of  the  McCuen  code  not 
treated  in  the  present  analysis  is  the  al-ility  to  specify  different 
thermal  conductivities  in  the  axial  and  radial  directions. 

The  comparison  between  the  presently  imported  code  and  those  discussed 
above  is  summarized  in  Table  1,  where  W2E 'ABLATE”  is  used  to  designate 
the  former.  The  full  name  of  the  code  reported  here  is  Two-dimensional 
Charring  Ablation  Program.  The  Hurwicz  code  is  omitted  from  the  table 
due  to  its  analog  rather  than  digital  natvre  and  its  non-applicability 
to  rocket  engine  thrust  chambers. 


Table  1- 


Comparison  of  the  Present  Program  With  Some  of  the  Previous  2D 
Conduction  and  Ablation  Codes 


TABLE 

Transient 

Conduction 

TIGER 

TAP 

yes 

MeCUBN 

2D- ABLATE 

(presently  reported 
code) 

yes 

yes 

yes 

t 

Charring 

no 

no 

no 

yes 

Surface 

Ablation 

no 

no 

yes 

yea 

More  than 
one  Material  ' 

yea 

yea 

no 

yes 

Temperature 

Dependent 

Properties 

no 

yes 

yes  . 

yes 

Finite  Dif¬ 
ference  Method 

explicit 

exp. 

exp.  \ 

Implicit  alternating 
direction  method  (ADM) 

Numerical 

Stability 

conditional 

cond. 

cond. 

unconditional 

Normal  Flux  at  no 

Curved  Boundaries 

no 

.  yes 

yes 

Anisotropic 

no 

no 

yes 

no 

Conductivity 


© 

MATHEMATICAL  MODEL 

PRELIMINARY  REMARKS 

A  two-dimensional  transient  mathematical  model  is  developed  belcw  des¬ 
cribing  the  coupled  conduction  and  charring  ablation  mechanisms  which 
are  induced  in  the  walls  of  a  rocket  engine  thrust  chamber  during  engine 
firing.  It  is  assumed  for  the  presentation  that  the  reader  is  familiar 
with  the  descriptive  motivation  and  discussion  contained  in  Ref.  2  ,  the 
final  report  of  the  first  phase  of  this  program.  Thus  the  discussion 
presented  here  will  be  limited  to  those  areas  in  which  the  Phase  2  model 
differs  from  the  Phase  1  one-dimensional  model.  These  include  the  extra 
O  dimension,  the  thermodynamic  simulation  of  the  chemical  reactions 

occurring  within  the  ablative  material,  the  emission  of  all  momentum 
effects,  and  the  degree  cf  assumed  knowledge  of  the  material  properties 
of  the  charring  ablative  material. 

The  general  approach  taken  is  similar  to  that  of  Ref.  2  in  that  internal 
charring  is  treated  in  depth.  Both  heat  flow  and  gas  mass-flow  are 
handled  in  two-dimensional  cylindrical  coordinates,  however,  using  the 
assumption  of  axial  symmetry  to  restrict  their  spatial  variation  to  the 
sxial  and  radial  directions. 
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In  this  section  we  will  describe  the  major  features  of  the  mathematical 
model*  Further  details,  such  as  the  basis  for  the  thermodynamic  simula- 
tion  of  chemical  reactions  in  the  energy  and  continuity  equations,  are 
given  in  separate  appendices  at  the  end  of  the  report* 

The  equations  developed  are  presumed  to  hold  in  regions  bounded  by  quite 
arbitrarily  curved  boundaries.  In  practice,  the  overall,  multi-region 
configurations  are  limited  to  those  typical  of  the  walls  of  rocket 
thrust  chambers.  Similarly  the  boundary  conditions  utilized  reflect  the 
heating  and  cooling  mechanisms  which  describe  the  thermal  environment  of 
the  chamber  walls  during  engine  operation.  Figure  1  is  a  sketch  of  a 
typical  thrust  chamber  wall  configuration  in  cross-section,  along  with 
examples  of  the  heating  conditions  encountered  (in  Fig.  1  the  receding 
hot  gas  boundary  is  indicated  as  well  as  the  normal  and  tAngqy^.-t^i 
directions  asswned  for  the  material  regions). 

EKEBGY,  comHirarr,  AHD  recession  equations 

I 

Using  cylindrical  coordinates  (axial  and  radial),  the  energy  and  contin¬ 
uity  equations,  respectively,  can  be  written  in  the  following  fora  for  a 
charring  wall  material  (see  Appendix  A  for  Justification): 
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Tj<pj.eal  Thrust  Chaaber  Configuration  ftyi^  Boi 
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where  i  ranges  over  the  chemical  reactions  talcing  place  within  the 
material.  For  a  non-charring  material,  in  'which  no  reactions  and  thus 
no  gas  mass  flow  occurs.  Eg..  (  2  )  does  not  apply  and  Eg.  ( 1 )  reduces  to 
the  pure  conduction  equation: 


Using  the  following  expression  for  the  density  P  of  the  charring 
material  (see  Appendix  B  for  the  derivation) 


+<  Kl-rc,S-rr,B)Fsc-Fa.cPBc’ 


.ve  may  substitute  in  Eqs  •  ( 1 )  end  (  2  )  to  obtain  the  following  more 
convenient  expressions: 
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r  dP 


'J\/rc,34rr(O-0 


dP 

sc 

dT 


„  fgc  ~  p  fjSscVTl  dT 
v^sc  dec  dT  ^Sdec  sc  dT  JjJ  dT 


Continuity 

dF 

p«F  -~EZ 
Vres  dT 


DF  Dp  Sp.  — . 

Fc,  s+rr,  s'1)  *#• +  Fdec  “#  +  Fsc  “^1 


Dg  . 

2sr+|^(yy 


(5) 


(6) 


Consideration  of  a  momentum  equation  is  avoided  by  the  use  of  tiro 
assumptions .  Strictly  speaking  K<3..  (  4)  la  an  expression  for  the 
density  only  of  the  porous  solid  material  in  the  char  region  and  not  of 
the  entire  region.  The  missing  term,  the  density  of  the  gas,  is  assumed 
to  vary  so  slightly  vith  tins  that  the  derivative  of  this  quantity  can 
be  neglected  in  the  continuity  equation.  Its  elimination  from  the  energy 
equation,  on  the  other  hand,  is  effected  algebraically  as  shown  in 
Appendix  A.  The  only  remaining  dependent  variables  (other  than  those 
parameters  with  fcncnra  temperature  dependence)  are  T,  G^,  and  G  • 


The  second  basis  upon  which  momentum  considerations  can  be  eliminated  is 
the  assumption  that  the  orientation  of  the  gas  mss  flowrate  in  the  char 


region  is  in  the  direction  of  the  temperature  gradient  vector,*  rather 
then  the  pressure  gradient  vector,  as  indicated  in  the  following  equation; 

-  $  /  I?  <*  •  «) 

In  the  energy  equation,  the  gas  cracking  reaction  Is  simulated  in  the 
convection  terms  by  rewriting  ^  and  ^  as  ^  and  |j“  there  H 

is  assumed  to  be  a  known  function  of  temperature  and  includes  the  effec¬ 
tive  heat  of  cracking  as  veil  as  the  integral  with  respect  to  temperature 
of  the  gas  specific  heat  at  constant  pressure. 

All  physical  am  chemical  material  properties  appearing  in  Eqs.  {  1 )  — 

(  6 ),  i«e.,  Fpy,  Fgc,  K»  PC,  and  H,  are  assumed  to  be  known 

functions  of  temperature. 


In  the  case  of  surface  erosion,  one  more  dependent  variable  of  the  analy- 

t  *4 

sis  must  be  accounted  for,  namely,  the  changing  radial  position, 
y  *  f(x,*0j  of  the  receding  hot  gas  boundary  (see  Hg.  l).  The  fol¬ 
lowing  differential  equation  relates  surface  recession  to  the  surface 
erosion  mechanisms: 


df 

df* 


(8) 


where  i  ranges  over  the  species  of  gas  or  liquid  generated  at  the  sur¬ 
face  and  dees  not  include  those  generated  within  the  material  and  ejected 

*33ie  direction  of  flow  may  be  further  restricted  as  discussed  on  page  65 
below. 
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at  the  surface.  The  values  of  G^  are  determined  through  a  temperature  ^ 
copied  analysis  of  the  erosive  mechanisms  occurring  at  the  surface. 

*7%iin  is  discussed  below  in  the  section  on  initial  and  boundary  eundi-  ' 
tlons.  In  Eq.«  {  8  ),  p  is  the  density  of  the  decomposing  vail  material. 
The  mimiR  sign  caa  the  ri^ht  side  of  Eq..  (  8  )  is  required  because  of  the 
negative  orientation  of  the  mass  flowrates,  G^. 

TTTTTTAT,  Aim  BGQOTttHy  CGSJDHTIG^ 

Tn  order  to  complete  the  specification,  of  the  model,  given  in  part  by 

i 

Eq.fi.  (5  )  -  (  8),  vhich  determines  the  tenperature  distribution  through¬ 
out  the  wall  nl r  and  the  mass  Xlcwrate  distribution  in  the  ablative 

material,  and  in  order  to  calculate  the  positions  of  the  receding  hot  gas 
boundary  due  to  surface  erosion,  it  is  necessary  to  specify  initial 
values  for  these  parameters  and  boundary  conditions  vhich  describe  the 

thermal  exchange  between  the  vail  materials  and  their  immediate 

*  • 

environment. 

For  purposes  of  the  present  analysis,  all  vail  materials  are  taten  ini¬ 
tially  to  be  at  seme  uniform  temperature.  After  firing  commences,  the 
vail  wyTtortais  are  assured  to  be  heated  at  the  hot  gas  boundary  by  con¬ 
vection,  but  are  permitted  to  radiate  to  the  surroundings.  The  other 
boundaries  may  be  either  insulated  or  radiating  surfaces. 


In  addition,  allowance  is  made  for  the' possibility  of  a  known  heat  flux 
at  certain  portions  of  the  boundary  due  to  heating  or  cooling  by  the 
outside  environment  (e.g.,  aerodynamic  heating  by  the  atmosphere  during 
re-entry),  finally,  perfect  thermal  contact  is  assumed  at  interfaces 
common  to  two  adjoining  materials,  which  implies  continuous  temperature 
and  heat  flux  distributions  across  the  interface.  Heat  is  transferred 
through  the  walls  by  conduction  only  until  some  point  within  the  inter¬ 
ior  of  the  charring  material  reaches  its  minimum  pyrolysis  temperature „ 
After  subsequent  generation  of  gas  and  formation  of  a  porous  char  within 
the  material,  convective  cooling  takes  place  due  to  the  flow  of  gas 
through  the  char.  In  addition,  heat  is  lost  or  gained  through  the 
process  of  gas  generation  itself.  Similarly,  all  the  heat  brought  to 
the  wall  at  the  hot  gas  boundary  is  conducted  away  from  the  surface  into 
the  interior  until  some  exposed  surface  point  of  one  of  the  wall  mate¬ 
rials  (or  the  char  residue  in  the  case  of  the  charring  material)  reaches 
its  minimum  decomposition  temperature .  Then  seme  of  the  heat  is  absorbed 
at  the  boundary  as  the  surface  erodes  (or  perhaps  released,  depending 
■upon  the  reaction  that  takes  place) .  In  addition,  the  surface  heating 
is  diminished  b :•  the  so-called  "blocking  effect"  of  the  gases  ejected 
into  the  hot  gas  stream,  including  those  generated  both  at  the  surface 
and  within  the  charring  mate'  '  1. 
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The  qualitative  discussion  given  above  of  the  driving  mechanisms  for 
the  Initial  and  boundary  value  problem  under  discussion  can  be  expressed 
more  fozmally  in  summary,  as  follows: 


Initial  Conditions 

-  — -  — -f—  -  ■  ■■■■■-■  . 


T(x,y,0)  »  Tq 

(9) 

Gyfoy,*)  “  0 

for  T(x,y,<r)  *  Tpy 

(10) 

f(x,T)  «  fQ(x) 

for  T(x,  f(x,«r),T)  *  Tdec  , 

(11) 

where  T  ,  T,  ,  and  T  are  known  constants  and  f  (x)  is  a  known 
o  dec7  py  o' 

function  of  axial  distance.  After  thfe  beginning  of  pyrolysis,  Eq.  (10) 
can  also  be  treated  as  a  boundary  condition  by  restricting  attention  to 
those  points  (x,y)  in  the  charring  material  for  which  T(x,y,  t)  a  T  . 

JtrJ 

Boundary  Condition 


If  n  is  taken  in  the  outward  normal  direction  at  a  boundary  (see  Fig. 
l),  we  may  write 


Kll0  heff(Taw"T)  +  2  <%)Gi  > 


(12) 


where 


r. 


heff  " 


h  +  h  ,  +  h 
conv  rad  env 


(15) 


0  at  insulated  boundary 


Equations  (12)  and  (1?)  comprise  a  convenient  formulation  of  all  the 
boundary  conditions  on  the  perimeter  of  the  multi-material  vail  being 
analyzed.  For  example,  Taw,  bconv>  and  the  species  simulation  (vhere 
the  range  of  i  is  the  same  as  in  Eq.  (8))  have  physical  significance 
only  at  the  hot  gas  surface  exposed  to  engine  firing.  At  other  bound¬ 
aries  they  either  are  zero  in  Eq.  (12)  or  cancel  out.  For  example,  in 
the  case  of  a  surface  at  temperature  T  radiating  to  an  environment  at 
temperature  T  .  ve  have  kCC(nv  *  0,  the  radiative  term  in  Eq.  (13) 
is  expressed  by 


hrad 


«£f(T 


,4 


en 


«> 


T  — T 
aw 


(14) 


and  (if  there  is  any  environmental  heating)  ve  have 


^env 

"env  “  T  -T 
av 


(15) 


Thus,  when  Eq.  (13)  is  substituted  in  (12),  hconv  and  TftW  do  not 

appear.  In  Eq.  (12),  T  is  a  known  function  of  x,  K  is  a  known 

function  of  T  f^r  each  material,  and  the  ( AHj)  are  known  constants 

for  each  material.  The  G.^  are  determined  by  a  method  developed  at 

Rocketdyne  (see  Ref.  9  and  Appendix  c)  which  combines  many  of  the 

features  of  the  Surface  Kinetics  Method  discussed  in  Ref.  2  with  a 

diffusion  approach.  In  Eq.  (13),  the  term  h _  includes  the 

conv 


22 


following  modifications  of  the  basic  heat  transfer  coefficient,  h 

conv' 

due  both  to  simulation  of  intermittent  engine  firing  with  soakback  and 
to  the  "blocking”  or  "blowing  effect"  of  the  gases  ejected  at  the  hot 
gas  surface: 


^conv  **  Fon^hconv^  “  BCp„  Y  Gt  O 


]h  Eq*  (l6)j  Foa(T)  is  a  known  step  function  which  takes  on  constant 
values  between  0  and  1  in  order  to  o  .emulate  "on"  time  during  representa¬ 
tive  time  periods  of  intermittent  steady  firing  or  high  frequency 

pulsing  with  soakback,  end  h  is  a  known  function  of  x.  The 

cony 

summation  in  Eq*  (16)  is  taken  over  all  gaseous  (only)  species  ejected 
at  thu  surface,  including  those  generated  internally  as  well  as  at  the 
surface*  In  the  blocking  term,  B,  C  ,  and  the  n.  are  known  constants 
and  the  Hj  are  known  functions  of  temperature*  For  the  gases  gener¬ 
ated  at  the  surface,  the  are  obtained  as  described  in  Appendix  C 
and  for  the  internally  generated  gas,  we  have 

& 

rt  /»2\  /.mV 


K  ♦  9. 
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Interface  Conditions 


At  an  interface  separating  material  regions  I  and  J,  the  continuity 
of  temperature  and  normal  heat  flux  conditions  are  expressed  in  the 
following  form: 


tr  _ 

-i$ir 


dT| 

3S|. 


(18} 


where  the  minus  sign  results  from  the  opposite  senses  of  the  outward 
normal  directions  from  materials  I  and  J. 


The  complete  initial  and  boundary  value  problem  that  has  been  solved 
numerically  and  coded  for  automatic  computation  is  given  by  Eqs.  (  5  )  ~ 
(18)  along  with  subsidiary  equations  from  Appendix  C  for  the  mass  flow- 
rates  of  the  gases  and  liquids  generated  at  the  hot  gas  surface. 
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NUMERICAL  SOLUTION  OF  ABLATION  AND  CONDUCTION  EQUATIONS 


© 


GENERAL  APPROACH 

The  tvo-dimens ional  transient  nonlinear  boundary  value  problem  of  abla¬ 
tion  and  heat  conduction  defined  by  Eqs*  (  5  )  -  (lO)  vas  solved  numer¬ 
ically  vising  finite  difference  techniques*  The  second  order  energy 
equation  vas  discretized  in  accordance  with  a  virtually  second  order 
generalization  of  the  unconditionally  stable  implicit  alternating- 
direction  method  of  Feaceman  and  Each  ford.  Generalization  and  modifi¬ 
cation  of  this  method  vere  required  for  several  reasons,  including: 

(l)  the  nonlinearity  of  the  equation  and  boundary  conditions;  (2)  addi¬ 
tional  gas  generation  and  convection  terms  in  the  equation  due  to 

charring;  (3)  coupling  of  the  energy  equation  with  the  continuity  and 
• 

recession  equations;  (4)  handling  of  more  than  one  material;  (3)  the 
treatment  of  normal  temperature  gradient  conditions  at  curved  bounda¬ 
ries;  and  (6)  the  changes  in  boundary  geometry  due  to  surface  erosion 
and  recession*  The  first  order  continuity  equation,  on  the  other  hand, 
vas  approximated  to  first  order  in  time  using  a  backward  difference  and 
to  second  order  in  distance*  Because  of  the  backward  time  difference, 
the  resulting  set  of  difference  equations  in  each  time  step  can  be 
solved  explicitly  by  proceeding  from  the  vicinity  of  the  char  front  to 
the  hot  gas  surface  with  a  four  point  formula. 


© 
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The  recession  equation,  i.e.,  £q,  (8),  was  discretized  using  a  forward 
difference  approximation  for  the  time  derivative*  Dus,  however,  to  the 
iterative  method  described  in  Appendix  C  for  evaluating  the  and 
new  values  of  the  radial  position,  y  •  f(x,  t),  of  the  hot  gas  hound** 
ary,  the  time  difference  approaches  a  second-order  central  difference 
approximation* 

■» 

The  use  of  essentially  central  or  backward  time  differences  in  -the 
finite  difference  version  of  Eqs •  (  5  )  -  (8)  should  guarantee  that  the 
solution  procedure  will  be  unconditionally  stable  with  respect  to  the 
time  step  procedure  used  in  the  interior  of  the  material  regions*  (The 
stability  of  the  distance  step  procedure,  on  the  other  hand,  used 
within  each  time  step  to  solve  the  continuity  equation,  has  not  .been 
established.  As  a  matter  of  fact,  a  homogeneous  linearization  of  the 
difference  equation  has  .been  shown  to  be  unstable.) 


As  a  result  of  the  high  values  of  the  heat  flux  encountered,  at  the  hot 
gas  boundary  during  periods  of  steady  firing,  limitations  are  imposed 
upon  the  size  of  the  time  steps  utilized*  In  general,  however,  these 
limitations  are  not  nearly  as  severe  as  that  required  for  an  explicit 
discretization  of  the  energy  equation,  and  are  not  necessary  at  all 
during  periods  of  sparse  intermittent  firing  or  during  soakback. 


26 


o 


Uncoupling  of  the  equations  is  achieved  in  each  time  step  by  predicting 


first  the  new  position  of  the  receding  boundary  (when  recession  occurs), 
then  the  two-dimensional  temperature  distribution,  and  finally  the  mass 
flowrates  In  the  (barring  material*  Only  one  such  calculation  cycle  is 
performed  in  each  time  step.  The  only  iterative  procedures  utilized 
are  local  in  character,  such  as  the  recession  calculation  (see  Appendix 
C)  and  a  special  routine  used  In  the  computer  program  for  computing 
values  of  the  effective  heat  transfer  coefficient  at  boundaries  subject 
to  strong  heat  fluxes* 

Spatially  the  discretization  is  achieved  through  imposition  of  a  mesh 
on  the  multi-iaaterlal  region  of  interest*  Mesh  points  are  located  at 
the  intersection  of  the  radial  and  axial  meBh  lines  with  each  other 
(regularly  spaced  points)  and  with  the  boundaries  and  interfaces 
(irregularly  spaced  points)  •  Axial  mesh  lines  must  be  equidistant,  but 
the  spacing  of  the  radial  mesh  lines  can  be  variable.  The  latter 
feature  permits  ccncsatsation  In  areas  of  possibly  strong  axial  heat 
flux  (such  as  the  nozzle  throat  region)  * 

In  simulating  normal  gradients  at  curved  boundaries,  exact  relation¬ 
ships  are  ecployed  In  the  usual  fashion  to  relate  the  normal  and  tan¬ 
gential  derivatives  to  the  axial  and  radial  derivatives  in  terms  of  the 
slope  of  the  curve*  These  derivatives  are  used  In  conjunction  with  the 
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physical  as  simp  t  ion.  of  no  local  flux  along  the  curve  and  second  order 
spatial  differencing  to  obtain  an  expression  for  the  normal  derivative. 
In  the  recession  procedure,  normal,  recession  Is  actually  predicted  and 
then  is  converted  to  the  equivalent  radial  recession,  again  in  terms  of 
the  slope  of  the  receding  boundary. 

After  the  iterative  prediction  of  the  new  radial  positions  of  the  hot 
gas  boundary  (i.e.,  the  nev  intersection  of  the  boundary  with  the  radial 
mesh  lines),  these  values  are  smoothed  to  conform  to  the  overall  con¬ 
verging  and  diverging  character  of  the  Initial  boundary  geometry.  Then 
first  derivatives  of  the  nev  curve  are  obtained  numerically  at  the 
radial  positions .  Finally,  interpolation  is  performed  to  obtain  both 
the  nev  axial  positions  and  the  derivatives  at  these  positions.  This 
provides  the  basis  for  the  calculations  in  the  remainder  of  the  time 
step. 

3h  the  time  step  procedure  utilized,  all  systems  of  equations  generated 

* 

are  at  vorst  tridiagonal  and  are  easily  solved  directly  or  by  well- 
known  recursion  formulas. 

Explicit  details  of  these  numerical  procedures  are  given  in  the  follow¬ 
ing  sections. 

& 

A  system  of  linear  equations  is  tridiagonal  if  its  nonzero  coefficients 
Appear  only  on  the  main  diagonal  and  its  immediately  adjacent  diagonals. 
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MESH  PROCEDURE  AHD  FINHE  DlFHBHEHCE  ISFIKITIOSS 

At  each  level  of  time  a  rectangular  mesh  or  network 

of  lines,  x  =  xi  and  y  =  y^,  where  x_^  =  x±~±  +  ^i-i 

y,  »  J4y,  is  Imposed  on  the  wall  materials  in  the  x-y  plane#  A 
0 

typical  example  is  shown  in  Fig.  2  .  Intersections,  (x^y^),  of  the 
mesh  lines  are  called  regular  mesh  points,  and  intersections  of  mesh 
lines  with  boundaries  or  interfaces  are  called  boundary  or  interface 
mesh  points.  Boundary  and  interface  points  are ^usutiliy  irregular.  A 
regular  mesh  point,  which  is  not  on  the  boundary  cf  the  material  region 
in  which  it  lies  or  on  its  interface  with  another  region,  is  further 
termed  an  interior  point. 

To  approximate  the  continuous  distributions  desired,  discrete  tempera¬ 
ture  and  Bass  flowrate  distributions  are  obtained  which  are  defined 
only  at  mash  points  and  at  discrete  times  t^.  These  approximations 
are  provided  by  employing  finite  differences  to  replace  the  time  and 
space  derivatives  of  Eq.s .  (  5  )  -  (18)  at  the  mesh  points  •  Three-level 
second  order  accurate  differences  are  used  for  the  space  derivatives, 
centered  at  interior  points  and  off-centered  at  boundary  and  interface 
points.  Two-level  differences  are  employed  for  the  time  derivatives, 
either  first  or  second  order  in  accuracy  depending  upon  the  weight 
given  to  the  time  levels  used  for  the  remaining  terms  of  the  equations 


Typical  Mash  Configuration  and  Categorization  of  Resulting  Mash  Points 


f 


±n  which  the  tiae  derivatives  appear.  Equations  {  }  )  -  (  7  )  are  dis¬ 
cretized  only  at  interior  points  and  Eqs.  (12)  and  (18)  only  at  boundary 
and  interface  points.  Equation  (8),  being  only  one-dimensional  in 
distance,  is  differenced  at  the  axial  locations 


For  convenience  of  notation  in  defining  centered  space  differences  at 
an  interior  point,  (x^,y^),  ve  assign  this  point  the  aubscript  0 
and  the  four  adjacent  points  the  subscripts  A,  B,  C,  and  D,  as  shown 
in  Figs,  2  and  3  »  where  the  lettered  points  need  not  be  interior 
points  (Fig,  2  ) ,  We  further  let  d^  be  the  distance  between  points 
0  and  A,  dg  the  distance  between  0  and  B,  etc.  (see  Fig.  3). 
Ihen  centered  finite  difference  approximations  to  the  first  and  second 
axial  derivatives  are  defined  for  use  in  Eqs.  (5)  -  (8)  as  follows, 
xihere  F(x,y,T) .  is  an  arbitrary  function  and  where  *  dA  +  dc 

8,11  V's'A.c 


5^0  T  “  6lP0,k  "  ^dCFA^'l'^aA"a(PIbJifS^S,i^,pA,C 

9  k 

01  1  S*V  •  2Cd0FA3k-DA,CP0)k^AFC,k^A,C 


(19) 

(20) 


Thst  defining  fczssalas  for  centered  redial  differences  are  analogous  to 
those  of  Eqs,  (l?)  aad  (20.)  and  thus  have  been  emitted.  Using  the  sscse 
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point  arrangement,  the  off-centered  (to  the  left)  first  axial  differ¬ 
ence  at  point  A  (which  i  leeded  when  A  is  either  a  boundary  or 
interface  point)  is  defined  by  the  following: 


O 


s|*  ^  “  8i">FA,k  - 

•  Hc-4>M<cv^iv3/v  •  <a: 

Similarly  we  obtain  the  off-centered  (to  the  right)  difference  at  C: 


6y;F_  .  =  (8  -d_B*)Fn  . 

X  \jf  X  X  v  X  x 


-  Ca^A,iI-DloFo,k+(ic-dc)pa,kVPA,c 


Both  left  anu  right  off-centered  temperature  differences  are  required 
at  an  interface  point  separating  materials  of  different  thermal  conduc¬ 
tivities  for,  in  such  a  case,  the  derivative  is  in  general  discontinu¬ 
ous  at  the  interface.  Downward  and  upward  off-centered  radial  differ¬ 
ences  are  defined  analogously  and  are  required  for  similar  reasons. 

All  the  expressions  defined  above  are  accurate  to  second  order  and  can 
be  derived  by  manipulation  of  Taylor  series,  truncated  after  the  quad¬ 
ratic  term,  expressed  at  each  of  the  points  A,  B,  C,  and  D  about  the 
point  0, 


The  time  differences  used  in  Eq.s  •  (5)  -  (8)  are  defined  differently, 
being  formulated  as  two-level  rather  than  three-level  expressions.  To 
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approximate  the  derivative  with  respect  to  time  at  any  mesh  point  0 
til 

in  the  (k+l)  time  step,  i.e.,  in  going  from  time  to  time 

we  use  the  following  definition  (where,  for  convenience,  we  drop 
the  subscript  on  At): 


Si 


*0,k»l  "  F0.k 
At 


k  <  k  <  k+l 


(23) 


Ihe  derivative  in  Eq.  {Zu)  is  taken  to  be  at  same  time  Tg  between  Tfe 
and  T  ,  When  k  is  k  or  k+l  exactly,  then  Eq,  (23)  is  a  for- 
ward  or  backward  difference,  ^ep.sctively,  %rxl  provides  only  first 
order  accuracy.  Only  when  §  «  k+i  dots  the  expression  give  full 
second  order  accuracy  >  The  value  of  Jc  is  djtenoined  in  a  given  dis¬ 
cretization  by  the  weighting  factors  assigned  to  the  spatial,  differ¬ 
ences  au3.  the  other  tto-varying  tersss  of  th*  equation  at  each  time 
levels  If#  for  caeaapl®,  to©  spatial  differences  were  assigned  weight 
ft  at  tiro®  Tj,  and  weight  3.  at,  1 1m  r.  j ,  then  k  would  assist® 
the  k+l*  If  via  term  a  dlsrretisfttlon  stable  when  errors  do  not 

%vm  during  th«  <&lsulAtloa  required  for  its  solution  and  call  it  uncon¬ 
ditionally  stabl©  when  it  is  stable  for  any  size  time  step  used,  then 
generally  taccnditloaal  stability  is  assured  only  when  £  i  k+i.  In 
the  next  section,  through  use  of  the  simple  linear  two-dimensional 
transient  heat  conduction  equation,  the  various  possibilities  for  £ 
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are  illustrated  and  a  stability  analysis  is  indicated •  Included  is  a 
comparison  of  sane  of  the  commonly  used  difference  methods  with  the 
alternating-direction  method,  vhich  for  the  linear  case,  effectively 
yields  a  value  for  £  of  k+&. 


DIFFERENCE  SCHEMES  FOR  SOLVING  THE  LINEAR  CONDUCTION  EQUATION 


As  indicated  above,  the  linear  two-dimensional  transient  equation  of 
heat  conduction,  which  follows,  is  a  convenient  vehicle  for  illustrating 
some  of  the  many  two-level  finite  difference  methods  for  solving  para¬ 
bolic  partial  differential  equations  such  as  the  energy  equation 
expressed  by  Eq.  (  5 ) : 


(24) 


Most  of  the  commonly  used  discretizations  of  Eq.  (zk)  can  be  repre¬ 
sented  by  the  following  two -parameter  family  of  finite  difference  equa¬ 
tions  (for  vhich  we  assume  a  constant  mesh  size  h  «*  Ax  =»  Ay  and,  as  a 
matter  of  convenience,  omit  spatial  subscripting): 


1  Tk+rTk 


5  J IF-*"  ^Tk+1+s5X1+(l-r)6^(l-s)6X 


(25) 


Hew  r,  8  are  weighting  parameters  which  lie  in  the  closed  interval 
[0,1]  and  whose  values  determine  the  level  of  the  time  difference, 
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i.a .,  the  value  of  £  (as  defined  in  the  last  section)*  In  general., 
ve  may  write 

k  »  k  +  (x*s)/2  •  (26) 

Thus,  for  example,  when  r  ■  s  *»  0,  the  most  familiar  form  of  Eq.  (25) 
is  obtained,  the  forward  difference  equation.  In  this  case  we  have 
£  «  k  and  only  first  order  accuracy  is  provided  by  the  time  differ¬ 
encing*  The  forward  difference  approximation  is  also  referred  to  as 
explicit  because  each  application  of  it  contains  only  one  unknown 
(temperature  at  time  T^+1)  and  each  member  of  the  system  of  equations 
generated  in  each  time  stop  can  be  solved  individually. 

If  we  designate  a  numerical  scheme  for  solving  a  system  of  finite 
difference  equations  as  stable  (in  time)  if  errors  do  not  grow  during 
the  calculation  and  unstable  if  they  do  grow,  it  can  be  proved  that  the 
explicit  scheme  (r  a  s  a  0)  is  stable  if  and  only  if  the  following 
inequality  holds: 

•/  (27) 

(A  stability  analysis  of  Eq.  (25)  for  all  values  of  r  and  s  is 

given  in  Appendix  D.)  Inequality  (27)  is  often  written  in  terms  of 

the  modulus  R  a  2£L 
hz 


as  follows* 


Because  this  condition  is  not  met  for  all  values  of  At,  we  say  then 
that  the  explicit  method  is  only  conditionally  stable.  This  restric¬ 
tion  is,  in  some  cases,  more  damaging  than  the  first  order  truncation 

error  incurred  in  using  a  forward  difference.  In  fact  for  large  values 

a 

of  -r,  the  scheme  may  not  be  useful  because  of  the  excessive  number 
n 

of  time  steps  required  despite  the  brevity  of  the  calculation  in  each 
time  step. 

By  setfirg  r  =  s  =  1  in  Eq.  (25),  we  obtain  the  implicit  backward 
difference  analog  of  Eq.  (24),  which,  although  stable  for  all  possible 
values  of  R  (i.e.,  unconditionally  stable),  is  only  first  order 
accurate  in  time,  £  taking  on  the  value  k+1.  In  this  case  simulta¬ 
neous  systems  of  difference  equations  are  generated  in  each  time  step 
containing  as  many  as  five  unknowns  each.  Such  systems  are  very  time 
consuming  to  solve,  usually  requiring  iterative  procedures  for  most 
efficient  treatment. 

For  r  «  s  «  &  and  thus  £  a  k+$,  Eq.  (25)  reduces  to  the  well-known 
Crank-Nicolson  method,  which  in  addition  to  being  unconditionally  stable 
is  second  order  accurate  in  time.  The  simultaneous  systems  of  differ¬ 
ence  equations  generated,  however,  as  in  the  case  of  backward  differ¬ 
ences,  require  time  consuming  iterative  methods  for  solution.  As  a 
matter  of  fact,  it  is  evident  that  as  long  as  r  and  s  are  both 


nonzero,  the  systems  of  equations  generated  will  be  of  the  same  form, 
requiring  iterative  solution.  If,  on  the  other  hand,  at  least  one  of 
the  parameters  is  0  in  each  time  step,  the  systems  generated  will  be 
tridiagonal  and  admit  of  simple  recursive  solution  without  iteration 
(the  form  of  the  general  tridiagonal  system  of  equations  is  given  below) • 
In  addition,  as  shown  in  Appendix  D,  it  is  necessary  for  unconditional 
stability  that  k  i  k+$.  Only  when  k  a  k+i  does  the  time  difference 
provide  second  order  accuracy. 

The  only  difference  schemes  which  meet  all  three  standards  (of  gener¬ 
ating  tridiagonal  systems  of  equations,  satisfying  a  necessary  criterion 
for  unconditional  stability,  and  having  second  order  local  accuracy  in 
the  time  difference)  are  those  for  which  r  »  0  and  s  »  1  or  r  »  1 
and  s  m  0.  If  either  of  these  assignments  of  values  is  used,  however, 
and  permitted  to  remain  fixed,  the  resulting  scheme  still  will  not  be 
unconditionally  stable;  although  the  assignment  satisfies  a  necessary 
condition  (k  *  k+J),  it  is  not  sufficient  to  provide  stability  for  all 
values  of  R.  However,  it  is  shown  in  Appendix  D  that,  by  choosing 
r  ■  1  and  s  «  0  in  odd  time  steps  and  r  »  0  and  s  «  1  in  even 
steps,  the  difference  scheme  is  unconditionally  stable.  By  specifying 
r  and  s  in  such  a  manner,  Eq.  (25)  reduces  to  the  alternating- 
direction  method  of  Peaceman  ind  Racbford .  The  generalization  of  this 
scheme  (given  below)  to  Eq.  (5)  retains  most  of  the  advantages 


described  here.  The  systems  of  equations  generated  are  tridiagonal  and 
the  value  of  k  is  k+fc  (before  linearization  in  which  coefficients 
are  taken  at  time  t^)  c  Although  a  full  stability  analysis  has  not  been 
performed  for  the  generalized  scheme,  it  has  usually  been  found  in  deal¬ 
ing  with  the  conduction  equation  that,  if  a  difference  procedure  is 
stable  in  the  linear  case,  it  will  be  stable  for  the  nonlinear  as  well. 
As  indicated  earlier  in  the  report,  however,  stability  of  a  difference 
analog  of  vhe  conduction  equation  does  not  necessarily  extend  in  general 
to  the  coupled  boundary  procedure.  It  has  been  shown  (Refs.  10  and  ll) 
that  for  the  linear  equation  in  a  bounded  region,  with  values  specified 
on  the  boundary,  the  alternating-direction  method  coupled  with  the  known 
distribution  on  the  boundary  is  unconditionally  stable.  In  the  case  of 
boundary  heat  flux  conditions,  however,  stability  does  not  automatically 
follow  and,  in  seme  cases,  restriction  of  the  size  of  AT  is  required, 
as  discussed  later  in  the  report* 

The  results  of  the  above  discussion  of  the  four  schemes  obtained  from 
Eq.  (25)  in  the  linear  case  are  summarized  in  Table  2  . 

COiPUTATIONAL  PROCEDURES 

In  each  time  step  of  the  numerical  solution  of  Eqs.  ( 5 )  -  (18),  the 
major  computational  procedures  are  performed  in  the  following  order: 


Table  2  •  Comparison  of  difference  schemes 
for  solving  the  linear  conduction  equation 


Case 

Difference 

Scheme 

Stability 

Criterion 

Value 
of  k 

Order  of 
Time 

Difference 

2>8a0 

Forward 

(explicit) 

0  £  R  a  *  i 

h2  ^ 

k 

iBt 

laQal 

Backward 

0  SR  <® 

k+1 

iet 

x 

Crank-aicolson 

0  SR  s» 

k+t 

z^ 

r»l,s«0 

and 

r*0, Sal 

Alternating-direction 

(Peaceman-Rachford) 

OSRSo 

k+i 

e”4 
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first,  n ev  surface  positions  are  predicted  (when  erosion  occurs)  at  the 
hot  gas  boundary;  second,  the  temperature  distribution  throughout  the 
vail  materials  is  obtained  at  the  nev  time  level;  and,  finally,  gas 
mass  flowrates  of  the  generated* gases  are  calculated  vithin  the  charring 
material  (after  charring  begins)  •  The  numerical  methods  employed  in  the 
program  for  each  of  these  steps  are  described  in  the  following 
paragraphs* 


Prediction  of  Haw  Surface  Positions 

First  an  iterative  calculation  of  the  values  of  the  mass  flowrates  of 
the  gases  and  liquids  generated  at  the  hot  gas  boundary  by  decomposition 
of  the  surface  wall  materials  is  made,  as  described  in  Appendix  C.  Then 
the  following  difference  analog  of  Eq*  ( 8  )  is  used  to  predict  the  new 
radial  positin  of  the  receding  surface  for  each  mesh  line  in 

the  time  step,  -* 

-  -  <*  w  /  **  •  (29) 

In  Eq«  (29),  the  svmoation  index  1  ranges  over  the  species  generated 
and  spatial  subscripting  is  omitted  for  convenience  of  notation*  The 
Iteration  described  in  Appendix  C  is  performed  in  such  a  manner  as  to 
approximate  the  average  value  of  the  0i  and  P  during  the  time  step. 


O 


<1 


Thu3,  the  subscript  k+J  is  justified  and  ve  see  that,  except  for  the 
time  level  at  which  the  space  derivative  of  f  is  obtained,  the  left 
side  of  Eq,  (29}  is  a  central  difference. 

After  prediction  of  the  new  positions  during  a  time  step,  the  new  values 
are  smoothed  by  a  reordering  process  to  insure  the  monotonicity  of  the 
converging  and  diverging  portions  of  the  thrust  chanfcer  as  initially 
given  in  the  program.  In  most  cases,  reordering  will  coca*  close  to 
preserving  the  volume  of  material  lost  and  is  such  fester  then  any  of 
the  more  sophisticated  methods,  such  as  curve-fitting.  Moreover, 
because  generally  the  raw  predictions  will  themselves  exhibit  monotoni¬ 
city,  reordering  will  cause  fever  unnecessary  changes  to  be  made  in  the 
predictions  than  will  most  curve-fitting  procedures . 

Derivatives  at  the  new  vail  positions  are  approximated  after  smoothing 
through  differencing  as  defined  by  5qs.  (19),  (21 ),  and  (22),  the  latter 
two  being  required  whenever  recession  takes  place  along  the  right  or 
left  extreme  radial  mesh  lines,  respectively.  Finally,  linear  interpo¬ 
lation  and  extrapolation  of  the  new  radial  positions,  and  derivatives  are 
used  to  obtain  the  corresponding  positions  and  derivatives  at  the  new 
intersections  of  the  boundary  with  the  axial  mesh  lines. 
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Finite  difference 


To  illustrate  the  linearization  of  the  generalized  Peace^san-Rachford 
alternating-direction  method  used  to  discretize  the  energy  equation,  we 
rewrite  Eq.o  (  5  )  together  with  an  acccopar’^ng  symbolic  scheme  to  indi¬ 
cate  the  time  levels  at  which  the  differences  have  been  taken.  (Eere 
k  and  k+2  represent  even  levels  and  k+1  an  odd  level): 


WrSr- 


b2  T  d2T 
3x'  dy2 


i  &£V  dxT/'ST'\2  fit?' 1  mrr  a*  - 

y  mx&j*  \5yj  J  “  m  A  d£^y  v 


ODD:  (k)  (k)  (k)(k+l)  (k)  (k)  (k) 


(k)  (k)  (k+l)  (k)  (k) 


EVEN:  (k  (k  (k  (k+l)  (k+2)  (k+2)  (k  (k+l)2  (k+l)x  (k  (k  (k+l)(k  (k 
+1)  +1)  +1)  +1)  (k+2)  +1)+1)  +l)+2) 


For  example,  in  either  an  odd  or  an  even  time  step,  the  symbol  (k+l) 

above  means  that  the  difference  analog  is  taken  at  time  level 

the  symbol  (k)x(k+l)  means  that  a  product  of  derivatives  has  been 

2 

differenced,  one  at  time  level  and  the  other  at  (k+l) 

means  both  differences  in  a  product  taken  at  tine  level  ^  and  k 
is  used  as  before  to  indicate  a  time  level  between  k  and  k+l  at 
which  a  time  difference  has  been  taken. 


The  difference  analogs  of  Eq..  (30)  at  an  interior  point  0  (see  Figs- 
2  and  3  ),  in.  odd  and  even  tima  steps,  as  indicated  in  the  symbolic 


«fc  % 

scheme,  take  the  following  forms,  respectively  (where  the  spatial  sub¬ 
script  0  is  omitted,  for  convenience  but  is  understood  to  apply  to  each 
parameter  appearing) s 

Odd  Step. 


(po)erc,k 

*  VyVtjtf  *¥yV  <Vj  Vk  * 


(31) 


Even  Step, 


(PC) 


T  ~T 
k+2  k+1 


ef f , k+1  At 


t+1  6xTk+l  x  k+XJ  *  t+1 


6  T. 


8A.i 


(32) 


In  Eqs.  (31)  and  (32);  (PC)eff  is  an  abbreviation  for  the  coefficient 


dT 


of  ^  in  Eq.  (5  )  and  is  discretized  as  follow: 


J 


(ecW  ■ 

<[t!sol;rc,8+rr,s-1+rdeo('W,k^raas(I’ao>j£.i 

+  <ide=Fac(Tmax,k)6Tmi.(Fdec)k^1 


(33) 
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where  (F.)  is  defined  for  each  gas  generation  reaction  1  by 

imax  1  k-& 


T 


(F,) 


max 


Fi^Tmax.  Ip  ”gl^Tmax,  k-1^ 

T  -T 
max,  k  max,  k-1 


(34) 


We  use  6m  (F.)  In  the  discretization  rather  than  the  derivative 
xmax  1  k-$ 

of  the  known  function  F. (T  „  .)  despite  the  introduction  of  a  third 

X  IUSOC;  K 

time  level,  the  r.dded  compi  ter  storage  requirements,  and  the  reduction 
to  first  order  local  accuracy.  This  is  to  ensure  obtaining  the  gas 
generation  effect,  even  with  a  time  lag,  in  the  event  of  large  Jumps 
in  temperature  for  which  intervals  of  large  values  of  F* (T  v)  would 
be  skipped  entirely.  This  is  especially  likely  when  using  an  implicit 
method  due  to  the  possibility  of  using  large  time  steps  which  would  not 
be  permitted  for  an  explicit  method  because  of  stability  considerations. 
Also,  we  use  r^mx  as  the  argument  of  Fi  because  the  gas  generation 
reactions  are  not  reversible  and,  in  case  of  a  drop  in  temperature  at 
some  point,  no  more  generation  can  occur  there  until  the  temperature  is 
raised  back  to  the  maximum  thus  far  attained. 


The  use  of  5^/6^  and  d^h^/byT^ .  instead  of  H'  (T^)  and 

H'tTj^)  in  Eqs.  (31)  and  (32)  is  similar  in  motivation  to  the  use  of 

(F.)  .  It  ensures  the  inclusion  of  all  heat  absorption  effects 

lwsx  1  k-| 

when  the  distance  increments  used  are  so  large  that  intervals  of  large 
H'(T)  might  be  skipped.  In  this  case,  moreover,  we  are  able  to  retain 
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second  order  local  accuracy  and  the  storage  requirements  are  not 
increased* 


It  is  evident  from  the  definition  of  the  difference  operators  &  and 

p 

6  that, t for  each  interior  point  0,  Eq.  (51 )  is  a  linear  algebraic 
equation  in  the  three  unknown  temperatures  Tq  and 

and  Eq.  (52)  is  linear  in  the  unknowns  TQ  k+2,  and  Tg 

Furthermore,  since  for  all  interior  points  0,  the  adjacent  points 
A,  D,  C,  and  D  will  ultimately  include  all  boundary  and  interface 
points  as  veil  as  the  interior  points,  we  will  be  left  with  many  more 
unknowns  tfcta  equations.  In  order  to  make  up  this  deficit,  we  obtain 
difference  analogs  of  the  boundary  and  interface  conditions  at  each  of 
the  respective  boundary  and  interface  points.  The  procedure  used  to  do 
this,  particularly  in  the  case  of  curved  surfaces,  is  the  subject  of 
the  next  section. 


As  indicated  above,  difference  analogs  are  required  of  the  normal 
gradient  conditions  expressed  by  EqB.  (12)  an3  (18)  at  the  boundary  and 
interface  points  in  order  to  generate  as  many  equations  as  there  are 
unknown  temperatures .  In  obtaining  the  added  equations,  we  must  be 
careful  not  to  introduce  any  further  unknowns  to  the  system,  or,  if  we 
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do,  to  provide  a  means  for  their  elimination*  In  both  Eq.  (12)  and 
Eq.  (is),  we  may  rewrite  the  normal  derivative  by  making  use  of  the 
following  exact  relationship,  obtained  through  geometrical  considera¬ 
tions  : 


„  .  &  3t  ,  /  *  r.  (tort 

3n  £  v3x  3x  cSyy  /  ^1+ 


(35) 


Here  y  a  f(x,r)  is  the  equation  of  the  boundary  or  interface  under 
consideration  (end  can  only  be  time  dependent  at  the  hot  gas  boundary), 
and  the  plus  (minus)  sign  is  used  at  a  lower  (upper)  boundary,  i.e,, 
when  the  outward  normal  has  a  positive  (negative)  component  in  tL* 
y-direction.  We  see  that  Eq*  05)  holds  even  when  the  slope  is  infinite 
by  taking  the  limit  as  ^  approaches  + 

At  r,  boundary  point  A  formed  by  the  intersection  of  a  horizontal  mash 
■  with  a  right  lower  boundary  (see  Pig.  4  ),  for  example,  we  substi¬ 
tute  i'Ac  05)  in  Eq.  (12),  obtaining 


(W  -**3  + 1  <®iW  -  Ki 


".a. -si; 


A  *  A 


(36) 


1  + 


<sj 


V  ^ 

To  discretize  Eq.  06),  we  |~j  to  second  order  using 

6i  'TA'  the  left  off- cents: red  $$  givpis  by  Eq.  (21 ).  We  have 
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no  similar  expression  for  4~f  ,  however  •  If,  tor  example,  we  express 
dg>| 

I  about  the  adjacent  intsric.v  point  0,  which  is  in  general  its 

'■V'A 

only  adjacent  interior  point,  we  obtain  (assuming  T  to  be  quadratic 
«**.  x)  •  *  * 


dTl  5t{ 


d2T 


I  lA  *  #l0  +  4a 


(57) 


We  eee  that  our  position  is  nos  improved  by  Eq.  (37)  due  to  the  appear¬ 
ance  of  the  mixed  derivative  at  0,  for  which  no  difference  equivalent 
has  been  defined  (the  reason  is  given  below) , 


The  difficulty  encountered  at  a  right  boundary  point  A  in  expressing 

3”j  or,  equivalently,  in  expressing  ,  can  be  circumvented  in 

several  ways  •  Those,  however,  which  provide  only  a  first  order 

Imatlon  to  or  are  not  compatible  with  the  alternating-direction 

method  (such  as  discretizing  4-3-1  directly)  were  not  considered*  A 

0x07,0 

second  order  procedure  for  evaluating  ,  patterned  after  a  method 
devised  by  R«  V*  Vlswanathan  for  solving  Poisson's  equation  by  relaxa¬ 
tion  (see  Ref*  12),  was  attempted  in  the  program  but  was  discarded  due 
to  the  poorly  conditioned  system  of  difference  equations  It  generated. 
This  method,  although  not  applicable  here,  might  be  quite  useful  when 
employed  in  conjunction  with  other  schemas  for  solving  the  energy 
equation*  A  complete  description  of  Its  use  in  connection  with  heat 
conduction  problems  can  be  found  in  Refs.  13  and  14. 
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The  method  finally  employed  to  simulate  normal  gradient  conditions  at 
curved  boundaries  requires,  In  addition  to  the  expression  given  by  Eq. 
(35)  for  ^  at  a  curved  surface,  a  similar  expression  for  ,  the 
counterclockwise  tangential  derivative  and  a  physical  assumption  about 
the  local  behavior  of  .  The  first  is  the  following  exact  egression 
relating  ^  to  the  axial  and  radial  derivatives: 


The  physically  based  assumption  is  the  vanishing  of  the  tangential 
derivative  at  a  curved  surface 


which  is  appropriate  for  most  situations  encountered  in  rocket  engine 
applications  due  to  the  ut.ua!  doaiiiiaticu  of  ^  by  *  Then,  from 
Eqs»  (38)  and  (39),  we  obtain  the  following  condition,  which  we  use  in 
2q»  (35)  to  eliminate  ^  in  an  odd  (even)  time  step  at  a  bound* 

ary  or  Interface  point  lying  on  a  horizontal  (^rii.eal)  mesh  line: 


dT  df  dT 


(40) 


In  an  odd  step,  we  combine  Eqs.  (35)  and  (40)  to  obtain 


and,  In  an  even  step,  ve  obtain 


where  the  signs  are  employed  as  In  Eq.  (35)# 

Using  Eqp#  (4l)  and  (42),  we  then  obtain  the  following  conditions,  equiva¬ 
lent  to  EqSo  (12)  and  (18),  for  use  at  boundary  and  Interface  points  on 
horizontal  mesh  lines  in  odd  time  steps  and  on  vertical  mesh  lines  In 
even  steps: 


Odd  Step#  At  a  right  or  left  boundary  point,  we  have 

+ 1  /\ 1  +  d§)2  •  (13) 

1  ’ 

/ 

At  an  Interface  point,  we  see  that  the  derivative  drops  out  entirely, 
yielding 


where  the  minus  sign  drops  out  o£  Eq#  { 18)  because  a  non-vertical  Inter¬ 
face  Is  both  a  lower  and  an  upper  boundary. 


Even  Step*  At  a  boundary  ve  have 


- K  |  -  ±  W  ♦ 1  /\| 1  + 


(45) 


and,  at  an  interface, 


h.  llj  -  kj  if. 


(46) 


dr 


We  see  that,  aside  from  the  physical  a&sisnption  that  ^  vanishes, 
Eqs.  (43)  -  (46)  are  exact* 


The  discretization  of  Eqs.  (43)  -  (46)  then  follows  quite  readily.  In 
an  odd  step,  -»  we  obtain  the  following. 

Left  Boundary  Point  C. 

“  i‘'Slc^hef^r  n~\k+ll 


C,£ 


Right  Boundary  Point  A 

,(-L  .  df 


V*  W  “  itlAjkl(WA>£C(WA“TA,^ 


+ ?  \ 


1  + 


<*u 


(47) 


1  +  ^lc,p2  ' 


(48) 
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Interface  Point  p  (with  region  I 


on  the  left  and  J  on  the  right) • 


(49) 


In  an  even  step,  -*•  t^,  the  following  equations  result: 

Lower  Boundary  Point  D. 


(+) 

"  «>W  K*J  -1B  ,3 

Dyk+1  D 


+  £(«,}  «g  }Ai  +  (Jfi  j 

i  D  2>^+1  'J  “*d,j «-r 


(  CT>\ 

\  /V/ 


tfeper  Boundary  Point  R. 


'  Wr  W  '  - 


♦  £<flg«g  3/ 

l  2  b#£+i  M 


i  ♦  (M  J  . 

'B^+r 


(51) 


Interface  Point 

(Kj) 


Q,3&fl  y  Q>k+2 


%^v2 


(57) 


ae  use  of  5  in  Eqsc  (47)  -  (52)  Indicates  that,  under  certain  heating 
conditions,  local  iterations  are  performed  in  order  to  obtain  averaged 
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values.  Except  for  the  case  of  surface  erosion  which  is  discussed  in 
Appendix  C,  these  iterations  are  performed  using  the  same  boundary 
formula,  but  explicitly  discretized.  A  new  value  of  temperature  is 
predicted  at  the  boundary  point  and  then  hQff  and  K  are  recomputed 
using  an  average  of  the  predicted  temperature  and  its  value  at  the  old 
time  level.  For  the  present  erslon  of  the  computer  program,  the  last 
tern  of  Eq.  (51)  vanishes,  surface  erosion  only  being  permitted  at  the 
hot  gas  boundary. 

In  an  odd  time  step,  then,  we  see  that  Eq.  (31)  with  Eqs.  (47)  -  (49) 
provides  complete  linear  systems  of  algebraic  equations  (one  for  each 
horizontal  mesh  line)  with  as  many  equations  as  unknowns.  Similarly 
Eqs.  (32),  (50)  -  (52)  generate  ccoiplete  systems  on  vertical  mesh  lines 
in  even  steps. 

Only  two  steps  in  the  method  used  to  obtain  values  for  temperature  at 
all  the  mesh  points  in  each  time  step  remain  to  be  described.  First, 
we  will  treat  the  trldiagonalizatlon  of  the  systems  of  equations  gener¬ 
ated  by  the  above  procedures  and  the  recursive  solution  of  each  trans¬ 
formed  system.  Secondly,  a  brief  discussion  is  included  of  the 
\ 

methods  used  to  calculate  the  temperatures  not  obtained  by  the  recur¬ 
sive  solution,  that  is,  the  temperatures  at  the  irregular  points  on 
vortical  (horizontal)  mesh  lines  in  odd  (even)  time  steps  (in  the 
case  illustrated  by  Fig.  2  ,  for  example,  we  require  additional  pro¬ 
cedures  to  obtain  in  even  steps) . 
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Tridiagonalizatlon  and  Recursive  Solution 

In  order  to  take  advantage  of  simple  recursive  procedures  for  solving 
systems  of  linear  equations,  we  wish  to  transform  the  systems  generated 
by  the  procedures  discussed  above  into  tridiagonal  systems  of  equations. 
A  tridiagonal  system  is  one  whose  coefficient  matrix  contains  nonsero 
elements  only  along  the  main  diagonal  and  its  immediately  adjacent 
diagonals.  Such  a  system,  say  consisting  of  n  equations  in  the  n 
unknowns,  u^.  i=l,...»n»  can  be  characterized  as  follows: 


'lul 

+  Clu2 

3  dl 

‘2*1 

+  b2u2  •+  c^ 

"  d2 

^2  +  b3u3  +  C3U4 

d3 

aiVi +  Yi +  biVi 

“dA 

a  ,u  „  +  b  ,u  , 

+  c  ,u  =»  d  , 

n-1  n-2  n-1  n-1 

n-1  n  n— 1 

Vn-l 

+  b  U  ad 

n  n  n 

The  solution  to  system  ( 53)  is  obtained  by  the  following  recursion  echeme 
(also  defined  in  Refs.  3  and  15)»  which  can  be  easily  derived  using  a 
direct  Gaussian  elimination  scheme  such  as  the  compact  Crout  method* 
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For  a  given  mesh  line  in,  say  an  odd  time  step  we  see  that  the  system 
of  equations  generated  by  Eqs.  (31)  and  (47)-(49)  complies  with  the 
conditions  for  tridiagonality  as  given  by  system  (53)  at  each  interior 
point  (using  Eq»-(3l))  and  does  not  at  either  boundary  point  (points  1 
and  n  of  system  (53))  or  at  an  interface  point.  If,  say,  the  ith 
point  on  the  mesh  line  is  an  interface  point,  then  the  system  generated 
above  has  the  following  form  (for  n  points)* 


o 
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To  eliminate  i^_2  and  i^+2  from  the  ith  equation  we  observe  that  for 
the  system  of  three  equations 

V  +  &i“  +  v  -  "1 


a2v  +  P2W  +  y2x  +  62y  +  *2Z  *  a2 

YjX  +  6^y  +  CjZ  =>  (r ^  » 

the  eliminant  to  remove  v  and  z  '  may  take  the  following  form* 

♦  aL(62e3-63e2)y 


(58) 


(59) 


Application  of  Bq.  (57)  twice  and  Eq.  (59)  for  as  many  interfaces  as 
appear  in  system  (55)  (only  one  is  shown)  will  transform  system  (55) 
into  a  tridiagonal  system  as  given  by  system  (55)*  Then  the  system  is 
solved  using  the  recursive  procedure  given  by  Eqs.  (54),  This  is  carried 
out  for  eaoh  successive  mesh  line,  horizontal  in  odd  time  steps  and 
vertical  in  even  steps. 


Obtaining  Temperatures  at  Irregular  Point w 

In  order  to  complete  the  temperature  calculation  in  eaoh  step,  a  sepa¬ 
rate  method  was  devised  and  coded  for  obtaining  temperatures  at  the 
irregular  boundary  and  interface  points  in  alternate  time  steps,  i.e„  at 
irregular  points  on  vertical  (horizontal)  mesh  lines  in  odd  (even)  time 


steps.  This  procedure  entails  the  use  of  linear  and  "constant"  inter¬ 
polation  and  extrapolation  (within  a  single  material  region  only)  to 
obtain  temperatures  at  the  irregular  points  from  those  already  calculated 
at  the  new  level  in  the  tridiagonal  systems  described  above.  At  all 
interface  points,  linear  or  constant  extrapolation  is  performed  from 
both  sides  and  the  average  is  taken.  Linear  extrapolation  is  used  where 
possible,  when  at  least  two  successive  interior  points  are  adjacent  on 
the  same  side  of  the  interface  point,  and  constant  extrapolation  if 
only  one  interior  point  is  available  (in  Pig.  5  ,  in  an  even  time  step, 
for  example,  linear  extrapolation  is  used  from  the  left  at  point  P  and 
constant  extrapolation  from  the  right).  Two  types  of  boundary  points 
are  considered:  those  exposed  to  engine  firing  and  those  not  exposed. 

In  the  latter  case,  extrapolation  from  the  interior  is  used  applying  the 
same  rule  as  at  an  interface,  but  without  taking  an  average.  At  an 
irregular  poiut  on  a  boundary  exposed  to  firing,  we  use  linear  surface 
interpolation  or  constant  surface  extrapolation  from  temperatures  already 
ealeulated  in  the  tridiagonal  systems  at  similarly  exposed  boundary 
points  of  the  same  material.  Thus,  taking  examples  from  Pig,  5  ,  in  an 

even  time  step,  Tg  is  extrapolated  from  Tg  and  Tg  interpolated 
5  4  2 

from  T-  and  T_  ,  while,  in  an  odd  time  step  T„  is  extrapolated 
1  3  B1 

from  Tg  and  both  Tg  and  Tg  are  interpolated  from  Tg  and 
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Because  the  method  adopted  for  irregular  points  does  not  directly 
utilize  any  of  the  equations  comprising  the  model ,  some  further  dis¬ 
cussion  of  the  approach  would  seem  to  be  in  order.  Apart  from  the  fact 
that  the  method  has  worked  in  the  cases  attempted,  it  should  be  noted 
that  it  does  simulate,  to  what  is  usually  a  first  order  approximation, 
the  thermal  behavior  of  adjacent  points,  which  was  in  turn  obtained  from 
the  model  equations.  Care,  however,  should  be  taken,  in  selecting  the 
mesh  lines  for  the  analysis,  to  prevent  extrapolation  over  relatively 
long  distances.  This  is  especially  true  in  the  uae  of  surface  extra- 
polaticn  at  boundaries  exposed  to  firing. 
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Discretization  and  Solution  of  the  Continuity  Equation 
Coupled  to  the  solution  of  the  energy  equation  (5)  for  the  temperature 
distribution  is  the  solution  of  the  continuity  equation  (6)  for  the  mass 
flowrate  of  the  generated  gases  in  the  charring  material.  After  sub¬ 
stitution  of  the  simplifying  assumption  given  by  Eq.  (7)  into  Eq.  (6), 
the  continuity  equation  takes  the  following  forma 


3x  +y  Jy  3t  * 


(60) 


where  p  is  the  density  of  the  charring  material  whose  time  derivative 
Is  given  by  the  left  side  of  Eq.  ( 6 ). 


O 
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In  the  numerical  solution  of  the  basic  equations,  the  mass  flowrate 

~ _ _  L 

analysis  is  uncoupled  from  the  temperature  calculation  bv  flolvlpg  tfrmT — - 

.  **• 

energy  equation  first  in  each  time  step  and  t^enTialng^he-caloulatsfi^  _  _  , 

4  ** 

values  of  temperature  in  the  solution  of  the'  continuity  equation*  Thirf4 
ie  achieved  in  the  discretization  of  Eq.  (60 )  by  expressing  the  time 

"  "  ■*  'V'.  ,  . 

derivative  as  a  backward  difference*  Spatially  the  discretisation  ~ 
somewhat  simplified  in  that  the  mass  flowrate  analyFta^Lguerf  ormsd 

only  at  regular  points.  Thus,  the  following  discretisation  of  Eq.  (| 

'  •» 

obtained  in  terms  of  mass  flowrates  and  temperatures  at  the  four  points 

6114  ^xi4-l»yj+i^  ^or»  aore  conveniently, 
the  pointe  1,  2,  3,  and  4  as  pictured  in  Fig.  6),  is  second  or$3r  in 
distance  and  first  order  in  time  and  ie  assumed  to  hold  at  the  point 

w+rttp  attl“  V 
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where,  for  convenience,  the  subscripts  y  and  i  have  been  dropped 
from  G  and  As.,  respectively.  In  Eq.  (6l),  for  each  reaction  i, 

y  * 

F.  ,  .  and  F.  e  .  ^  are  evaluated  at  T  .  .  and  T  ,  . 

i  f  X  f  k  i)X|  max  i  f  k  id&x  f  X  f  Ji-*! 

respectively,  as  in  Eqs.  (33)  and  (34).  For  each  set  of  points  1,  2,  3, 
and  4,  such  that  at  least  point  2  is  in  the  charring  material,  Eq.  (61 ) 
is  rewritten  in  the  following  form  and  solved  explicitly  for  G^  ^ i 

23+1  "  R2,k  A*^  G2,k  =  ^  +  2j+l  *  R4,k  A^G4,k 

+  ^  +  23+1  “  'S.k  ^C3,k"^  “  23+1  +  Rl,k  Ax^Gl,k  (62) 

4 

j&r  .(*)- 
“2  Z  6t  pl,k 
1 


where  ^  is  a  symbolic  representation  of  the  term  of  the 

summation  in  Eq.  (6l).  Equation  (62)  is  used  to  sweep  across  each 
horizontal  mesh  line  from  left  to  right  starting  with  the  topmost,  mesh 
line  and  going  down  from  line  to  line.  The  boundary  conditions  are 
simulated  by  the  following  interpretations  of  Eq.  (62) 2 
a)  If  an  interface  common  to  the  charring  material  and  one  of  the  other 
materials  cuts  through  the  rectangle  formed  by  points  1,  2,  3,  and  4, 
then  both  G  and  F^  vanish  at  all  points  in  the  other  material.  If 
®ny  of  the  points  l  (other  than  point  2)  happens  to  be  on  the  interface, 
then  G.  is  set  equal  to  0,  but  ?<  .  can  still  be  nonzero.  If  point  2 

x  1#* 
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is  on  an  interface  (in  which  case  at  least  one  of  the  other  points  must 

he  in  the  charring  material)  „  both  G  g  and  g  can  he  nonzei . . 

b)  If  a  boundary  cuts  through  the  rectangle,  again  both  G  and  F^ 

are  zoro  at  points  outside  the  material,  and  again  F^  can  be  nonzero 

on  the  boundary.  G  is  zero  at  an  upper  boundary,  but  can  be  nonzero 

at  a  nonius  ulated  lower  boundary,  c)  If  either  T _ „  or  T  .is 

mai,^  max, 4 

less  than  T  or  if  point  4  is  outside  the  charring  material,  then  G0 
ET  * 

is  set  equal  to  0.  If  G.  is  0.  then  G0  is  obtained  from  the  following 
two-point  formula  instead  of  £q.  (62)2 


2,k 


5m(.  +  % 

i  4,k  Ai  vp2,k*  p4,k^ 


(63) 


In  practice,  a  .modification  of  £q.  (?)  is  used  to  obtain  2^  k  at  each 

point  l  for  use  in  £q.  (62).  Because  we  know  that  the  direction  of  the 

mass  flowrate  -rector  to  the  right  of  the  throat  plane  should  be  downward 

to  the  right,  we  limit  the  range  of  allowable  values  of  ^  to  the 

closed  internal  I— 1 ,0].  ^  is  obtained  at  first  as  the  ratio  of 

6jT,  .  to  fi  D,  ..  and  then  is  (a)  set  to  zero  if  the  ratio  is  positive, 

%  x,&  y  z,& 

(h)  set  to  -1  if  the  ratio  is  less  than  -1,  and  (c)  left  unchanged  other¬ 
wise.  Do  the  left  cf  the  throat,  we  similarly  restrict  2^  k  to  the 
rarre  [0,l]  with  one  exception.  Because  2^  ^  is  nonnegative  to  the  left 
of  the  throat  we  must  be  sure  that  the  coefficient  of  Gg  ^  in  £q.  (62) 
does  not  vanish.  This  is  done  in  the  computer  program  by  automatically 
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calculating  the  minimus  value  ^ 


O 


1 

=  (l-  2j^')  for  all  values  of  J 
and  fr.  to  the  left  of  the  throat  and  then  restricting  R.  .  to  the 
range  [l,r]  where  r  =  min  {f^,!},  where  0  <  f  <  1. 

The  temporal  stability  of  the  mass  flowrate  solution  procedure  given  by 
Eq.  (62)  is  assured  due  to  the  backward  time  difference  employed,  provided 
that  tho  schemes  used  to  solve  both  the  energy  equation  and  its  boundary 
conditions  are  also  stable.  The  stability  with  respect  to  distance  of 
the  procedure  as  used  explicitly  in  each  time  step  from  point  to  point 
is,  however,  another  matter.  A  complete  stability  dialysis  of  the  latter 
type  has  not  been  performed,  but  it  can  be  shown  (by  methods  similar  to 
those  used  in  Appendiz  0)  that  a  homogeneous  linearization  of  Eq.  (62) 
is  unconditionally  unstable  in  the  y  direction  and  in  the  x  direction 
to  the  right  of  the  throat.  To  the  left  of  the  throat,  on  the  other 
hand,  it  is  unconditionally  stable.  These  results  should  not,  however, 
bt  surprising  because  we  know  (due  to  the  restrictions  on  R(x,y,T))  that 
the  solution  must  grow,  to  the  right  of  the  throat,  in  the  same  axial 
and  '  -  ’■’nl  directions  in  which  the  numerical  solution  is  marched  out; 
whez»,  .  the  left  of  the  throat,  the  solution  must  grow  in  the  same 
radial  direction  but  in  the  opposite  axial  diroction.  No  steps  have 
been  taken  to  introduce  at  least  conditional  stability  to  the  numerical 
procedure  because  in  practice  it  has  proved  to  bo  quite  satisfactory 
when  used  in  conjunction  with  the  temperature  calculation.  Only  a 
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\ 

relatively'^  f  ew  distance  steps  are  taken  in  performing  the  calculation 

and  the  varies  generated  are  probably  small  enough  that  local  instability 
t 

will  not  affect  the  overall  stable  time  step  calculation.  In  any  oase, 
it  is  felt  ty  the  writers  that  no  usable  modification  of  the  numerical 
scheme  would'.provide  stability  in  the  axial  direction  to  the  right  of 
the  throat  du^  to  the  sign  of  R(x,y,?),  now  would  any  be  stable  radially 
on  either  sid^  of  the  throat  because  we  must  march  from  higher  to  lower 

i 

I 

values  of  j  jlue  to  the  available  boundary  conditions. 

♦ 

» 

; 

Control  of  Size! of  Time  Inoremente 
♦ 

In  the  one-dimexteional  code  completed  in  .'Phase  1 .  of  the  present  program 
(see  Ref .  2),  th|  approach  taken  to  the  maximization  of  the  time  step 
was  to  start  with^  a  very  email  increment  and  allow  it  to  increase  period- 
ically  during  the  program  run  provided  that  several  built-in  criteria  are 

j 

satisfied'.  Two  options  are  available  in  the  two-dimensional  program,  one 

\ 

similar  in  approach^  to  the  Phase  1  program  and  the  other  rather  different. 

i  ' 

•  \ 

The  first,  vhioh  we .shall  call  Option  A,  consists  simply’ of  scheduling  a 
variable  array  of  tiiie-ateps  directly  ae  input.  (Option  A  was  exercised 
for  the  second  checkout  case  described  in  the  Results  Section  below  and 
is  given  in  Table  4.  ’Option  B,  discussed  below,  wes  used  for  the  firot 


checkout  case.) 


In  Option  B,  on  the  othj&r  hand,  a  single  largo  (input)  time  step  is 
employed  wherever  possible  throughout  the  program  run  and  is  only  reduced 
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temporarily  as  required.  Such  reductions  are  aoheduled  within  the  pro¬ 
gram  in  accordance  with  a  stability  eriterion  based  upon  the  intensity 
of  engine  firing  as  simulated  in  representative  periods  by  the  fraction 
of  "on"  time,  ?oq(t)  (defined  below  Eq.  (l6)). 

It  has  been  shown  by  previous  investigators  (see  Ref.  (16),  p.  47-49, 
for  example)  that  instability  can  be  introduced  in  a  numerical  solution 
by  the  intensity  of  heating  at  the  boundary.  For  the  one  and  two 
dimensional  forward  difference  methods  used  to  solve  the  heat  conduction 
equation,  the  following  inequalities  must  be  satisfied,  respectively,  to 
ensure  a  stable  solutionis 


R*^a<2Tu5T 


(64) 


E  <  2C24ST  ^ 

% 

where  N  3  (hAx)/K,  A*  =  Ay,  and  h  and  K  are  constants.  In  the 
present  program  the  following  analogous  criterion  is  used  to  limit  the 
size  of  At  when  Option  B  is  exercised: 


R  * 


ofrr  .  Q 
(to)2  **. 


where  As  «  min[oin(Ax.), 

i  # 

determined  constant.  In 


(66) 

Ar]»  «  (heffAz)/K,  and  C  is  an  empirically 
order  to  be  able,  within  the  program,  to  schedule 


*C  wan  taken  to  be  .4  or  .6  in  all  checkout  runs  using  Option  B. 
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•11  time-step  reductions  based  on  inequality  (66)  before  the  numerical 
calculation  begins,  a  reference  temperature,  equal  to  60#  of  the  throat 
value  of  T4w(x),  1b  used  to  compute  values  for  hoff  and  K. 

There  is  one  other  criterion  for  temporary  time-step  reduction,  based 
upon  amount  of  surface  recession,  which  is  employed  with  either  Option 
A  or  B  when  surface  erosion  takes  place .  The  purpose  is  to  prevent 
excessive  radial  recession  in  any  one  time  step.  This  is  accomplished 
in  the  program  by  specifying  as  input  a  fraction  of  the  radial  increment 
Ay  as  the  maximum  allowable  amount  of  radial  surface  recession  along 
any  mesh  line  x  =  x^  in  a  time  step.  When  this  criterion  is  violated, 
the  program  automatically  reduces  the  time  step,  through  multiplication 
by  an  input  constant  less  than  one  in  magnitude,  and  then  begins  the 
time  step  again  ( thia~does~n6T~iirTOlve-  much- recalculation  because  re¬ 
cession  is  the  first  mechanism  in  each  time  step). 

One  other  means  of  time  step  control  is  provided  in  the  program  by 
permitting  the  optional  partition  of  the 'first  time  step  (only)  into  aa 
many  smaller  steps  ae  desired.  This  may  be  dons  is  conjunction  with 
either  Option  A  or  B. 

An  excellent  feature  of  the  Phase  1  code  is  its  capability  of  varying 
the  distance  increment  with  time.  This  was  not  considered  feasible  in 
the  Phase  2  code  because  of  the  exoessive  amount  of  recalculation  which 


would  be  required  to  establish  the  new  mesh.  However,  in  the  two-dimen¬ 
sional  program,  added  flexibility  is  provided  by  the  variable  spacing 
permitted  on  the  radial  mesh  lines. 


DESCRIPTION  OF  THE  COMPUTER  PROGRAM 


GENERAL  REMARKS 

She  numerical  procedures  developed  'above  for  solving  the  two-dimensional 
ablation  and  heat  conduction  model  given  by  Fqs.  (5  )  -  (18)  were  coded 
in  Fortran  IV  for  automatic  computation  on  the  EDFM  7^9^*  The  computer 
program  has  been  given  the  designation  "2D-ABLATEs” 

Input  to  the.  program  is  handled  by  a  very  sophisticated  input  routine 
vhich  was  adapted  from  a  current  Rocketdyne  procedure*  Its  advantages 
Include  reading  all  input  into  a  single  data  array  (vhich  is  allocated 
within  the  program),  vhich  can  thus  be  conveniently  augmented  as 
required,  and  a  particularly  efficient  method  of  handling  multiple  runs 
or  updating  data  from  run  to  run  through  use  of  single  cards  to  replace 
the  contents  of  specific  locations  in  the  data  array. 

In  order  to  take  advantage  of  the  overlay  feature  of  the  Fortran  IV 
system,  the  program  was  written  in  several  links  •  In  addition,  to  the 
main  control  link,  Link  0,  two  other  links  are  provided  for:  (a)  Link 
1,  which  hardies  all  execution  prior  to  the  time-step  calculation,  and 
(b)  Link  2,  the  time-step  calculation,  vhich  carries  out  the  numerical 
procedures  described  earlier. 


Link  1  Includes  allocation  of  the  input  data  array,  imposition  of  mesh 
lines  on  the  input  curved  boundaries  and  Interfaces  resulting  in  the 
regular  and  irregular  .nesh  points  defined  earlier,  generation  of  des¬ 
criptive  arrays  characterizing  these  points,  linear  Interpolation  at  the 
positions  xi  m  x^_1  +  Axi  of  the  input  axially  varying  tabular  data, 
piecewise  linear  and  quadratic  fitting  of  the  input  temperature  dependent 
material  properties,  r  lecevise  constant  and  linear  fitting  of  the  time 
varying  data  (e.g*,  ?oa(r),  which  is  defined  earlier),  and  scheduling 
temporary  reductions  of  the  input  computational  time-step  (as  required 
by  Inequality  (66))  for  use  in  periods  of  strong  surface  heating  unless 
the  option  is  exercised  of  scheduling  variable  time  steps  directly  as  input. 

In  order  to  acquaint  the  reader  with  the  scope  of  the  program,  a  dis¬ 
cussion  follows  of  the  input  data  required  for  its  operation,  brief 
descriptions  are  given  of  the  execution  performed  in  Links  1  and  2, 
overall  flowcharts  are  presented  of  the  main  program  logic,  and  the 
output  options  available  in  the  program  are  outlined. 

INPUT  REQUIRED  FOR  THE  COMPUTER  PROGRAM 

The  input  data  required  for  operation  of  the  2D-ABLATE  computer  code 
include  the  following: 

(a)  Thrust  chamber  wall  geometry. 

Piecewise  quadratic  functions  are  used  to  describe  the  configuration  of 
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all  material  boundaries  and  interfaces*  Because  the  materials  are 
numbered  from  1  to  a  maximum  of  5,  each  quadratic  piece  must  also  be 
identified  by  the  number  assigned  to  the  material  or  materials  (in  case 
of  an  interface)  which  it  bounds*  The  kind  of  boundary  condition 
assumed  must  also  be  read  in;  possibilities  include  insulation,  engine 
firing,  radiation  to  the  surroundings,  environmental  heating,  or  sane 

combination  of  these*  A  maximum  of  kO  quadratic  pieces  may  be  used  in 

) 

describing  the  wall  geometry* 

(b)  Discretization  parameters. 

The  only  required  parameters  for  the  discretization  are  the  basic  time 
increment  At  (Option  B)  or  an  array  of  time  increments  (Option  A) 
scheduled  as  described  in  (g)  below,  the  radial  distance  increment  Ay, 
and  the  axial  increments  Ax±,  where  i  goes  from  1  to  a  maximum  of  15* 
It  is  not  necessary  to  read  in  the  coordinates  of  any  mesh  points;  these 
are  all  automatically  determined  in  Link  1*  Due  to  storage  requirements, 
£y  must  be  chosen  so  that  no  more  than  55  horizontal  mesh  Ixuss  will 
Intersect  the  wall  boundaries. 

(c)  Temperature  dependent  material  properties. 

The  product,  PC,  of  density  and  specific  heat  is  read  in  as  a  piece- 
vise  linear  function  of  temperature  for  each  wall  material  and  the 
conductivity,  K,  as  a  piecewise  quadratic.  The  effective  enthalpy  of 
the  gases  generated  internally  in  the  charring  material  is  also  read  in 
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as  a  piecewise  linear  function  of  temperature .  Each  of  the  re¬ 
generation  reactions  which  take  place  within  the  charring  ablator  (a 
maximum  of  three  is  permitted)  is  characterized  by  the  temperature  range 
in  which  it  occurs  and  by  the  maximal  value  of  F^(t)  (see  Honenclature 
or  Appendix  B)  which  is  attained  during  the  reaction.  These  values  are 
then  used  within  the  program  to  generate  a  trigonometric  function  which 
is  employed  to  describe  the  variation  of  F^T)  within  the  reaction 
ranges 

(d)  Constants  for  gas  generation  reactions  In  the  charring  material. 

The  constants  required  to  characterize  the  reactions  are  Fres, 

Ssc*  Slec*  Pv*  ps*  rc,s'  an*  ^r,s* 

i 

(e)  Hot  wall  erosion  parameters. 

Included  in  the  input'  data  array  is  a  list  of  all  the  parameters  required 

« 

to  perform  the  calculations  given  in  Appendix  C.  They  include 

(1)  flew  constants  of  the  CGehustloa  gas  2  specific  heat  ratio, 
molecular  weight.  Chamber  pressure,  throat  radius,  radius  of 
curvature  of  the  throat  boundary,  and  parameters  used  in 
computing  the  blocking  effect. 

(2)  melting  rate  constants  for  each  exposed  materials  thermal 
conductivity  of  the  melt,  viscosity  parameters,  specific  heat 
of  the  melt,  heat  of  fusion,  and  melting  teeperature. 


(5)  vaporization  constants  for  each  material :  moss  diffusion 
coefficient  at  standard  conditions,  molecular  veight  of  the 
vapors,  vapor  pressure  parameters,  heat  of  vaporization,  and 
density  of  the  material* 

(4)  chemical  reaction  constants  for  each  material:  molecular 
veight  of  reactant  gas,  mol  fraction  of  reacting  gas,  rate 
equation  constants,  enthalpy  change  for  the  reaction,  mass 
diffusion  coefficient  at  standard  coalitions,  and  specific 
heat  of  the  reaction  products . 

(f)  Parameters  affecting  the  heat  flux* 

Values  are  required  for  each  of  the  following  parameters  at  various 

prespecified  axial  positions:  T^,  h„oav,  0f,  and  q^.  In  addition, 

two  radiation  constants  are  required:  cr«  and  T _ . 

env 

(g)  Time  varying  parameters* 

To  characterize  intermittent  firing  and  pulsing  with  soakback,  an  array 
of  values  of  ?on(r)  la  read  in,  each  of  which  is  held  constant  during 
a  preassigned  time  interval.  If  Option  A  is  exercised,  an  input  array 
of  time  steps  is  similarly  associated  with  these  intervals*  An  array  of 
weighting  constants  is  also  read  in  to  permit  piecewise  linear  variation 
of  with  time.  Hone  of  these  arrays  may  contain  more  than  40 

elements* 

(h)  Initial  temperature* 


The  initial  wall  temperature  is  required  os  an  input  constant 


LIHK  i  PROCEOJRES 


Input  Processing  and  General  Procedure 

All  data  for  the  program  are  arranged  In  a  matrix  format  -which  is  input 
as  a  single  array  and  stored  on  tape  by  the  program  so  as  to  be  made 
available  if  needed  at  the  completion  of  the  problem  run*  Thus  to 
perform  multiple  runs  it  is  only  necessary  to  read  In  those  data  ele¬ 
ments  (along  vith  their  array  positions)  whose  values  change  from  run 
to  run* 

Error  checks 'are  made  extensively  in  Link  1  and,  in  the  event  of  an 
error,  execution  continues  In  an  attempt  to  check  for  other  errors 
prior  to  termination.  Data  Input  which  defines  the  geometrical  config¬ 
uration  is  especially  critical,  and  it  is  Important  that  certain  rules 
(discussed  in  the  next  section)  be  strictly  followed  in  choosing  axial 
and  radial  increments  to  be  used  by  the  program  In  imposing  a  mesh  on 
-the  wall  configuration.  The  mesh  Imposition  essentially  consists  of 
determining  intersections  of  the  functions  describing  the  input  bound¬ 
ary  and  interface  geometry  vith  the  horizontal  and  vertical  mesh  lines. 
Arrays  are  generated  and  stored  consisting  of  these  locations  along 
with  the  derivatives  of  the  boundary  and  Interface  functions  at  these 

i 

locations,  boundary  types,  and  region  numbers  associated  vith  them. 


o 


Temperature  dependent  data  are  than  processed,  Including  enthalpy,  PC, 
and  thermal  conductivity*  The  first  tvr>  sets  are  fitted  by  several 
linear  se^centa  over  the  temperature  range  -while  thermal  ccaaductivlty 
Is  fitted  by  quadratic  segments* 

Further  processing  Is  required  In  Link  1  to  linearly  lntezpola.tr  the 
axially  varying  Input  parameters,  T^  ^com3  and  ^  at  the 
axial  positions  of  -the  vertical  mesh  lines  • 

Then, -  unless  Option  A  is  exercised  of  scheduling  variable  time  steps  as 
Input,  &  schedule  of  computational  time  Increments  Is  generated  In 
accordance  vlth  the  stability  criterion  given  by  Inequality  (66)  above, 
"'which  depends  tpon  the  Input' values  of  *’on(<r)  used  to  simulate 
Intensity  of  engine  firing  In  prespecified  time  Intervals. 

* 

Finally,  Initialization  of  certain  Link  2  parameters  Is  performed*  This 
terminates  Link  1  and  program  control  Is  then  transferred  to  Link  2* 

A  flowchart  of  the  overall  Link  1  logic  Is  shown  In  Fig.  7  • 

i 

i 

Some  Buies  Governing  Thrust-  Chamber  Geometry  Definition 

1.  A  mesh  line  •Khlch  Intersects  the  configuration  must  contain  exactly 
-two  boundary  points  which  are  terminal  points  of  the  line* 

2*  A  mesh  line  may  contain  as  mazy  as  eight  boundary  and  Interface 
points* 


77 


vmm 


'(  Link  1  Procedure 
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(zero  out  if  new  prob¬ 
lem  or  Input  from  tape 
if  multiple  run) 


, - - 1 - 
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of  axial  varying 
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puted  at  verti¬ 
cal  mesh  lines 

fr:  v  v 
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mesh  upon  wall 
'  geometry  *■  d 
generation  of 
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terizing  regular 
and  irregular 
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generated  arrays 
on  signal 
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distance  varying 
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Generate  array 
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time  Increments 

based  on  F  (t) 
on' 

Initialize  vari¬ 
ous  arrays  and 
parameters  used 
in  Link  2 


End  of 
Link  1 


3 «  A  horizontal  mesh  line  may  contain  as  many  as  23  points  of  all 
kinds  o'  A  vertical  mesh  line  my  contain  as  many  as  43  points  0 

V  * 

X 

4.  A  mesh  line  which  intersects  the  configuration  must  contain  a 
minimum  of  four  points;  these  my  be  any  combination  of  boundary,  inter¬ 
face,  and  regular  mssh  points o 

5 «  At  least  one  regular  mesh  point  must  lie  between  successive  bound¬ 
ary  and  interface  points# 

6  o  No  common  intersection  of  three  or  more  quadratic  boundary  or 
interface  sogmc  its  my  lie  on  a  mesh  line# 

7.  No  mesh  line  my  coincide  with  a  vertical  or  horizontal  boundary 
or  interface* 

LINK  2  EXECUTION 

In  Link  2  of  the  program,  the  entire  time-step  calculation  performed 

* 

as  described  earlier  in  the  section  on  numerical  procedures.  The  pro¬ 
gram  flow  is  iterative  in  nature,  calculations  being  performed  in  one 
time  step  at  a  time  in  the  following  order:  an  iterative  procedure 
(see  Appendix  C)  is  used  to  determine  the  erosion  characteristics  at 
the  hot  ©as  boundary  and  predict  the  resulting  surface  recession,  the 
temperature  distribution  is  obtained  using  the  alternating-direction 


<4W. 


method,  and  the  gas  mass  flowrate  calculation  is  performed  for  gases 
generated  in  the  charring  ablator. 

Of  these  calculations,  the  temperature  computation  requires  the  most 
intricate  program  flow  due  to  the  need  for  alternately  converting: 
horizontal  arrays  (at  both  regular  and  irregular  points)  and  irregu¬ 
larly  located  vertical  arrays  at  odd  time  levels  to  and  from  the  cor¬ 
responding  vertical  and  irregularly  located  horizontal  arrays  at  even 
levels . 

A  flowchart  of  the  overall  calculation  is  given  in  Figure  8. 

OUTPUT  PROCEDURES 

Output  from  uhe  Phase  2  code  2D-ABLATE  is  printed  out  in  both  Link  1 
and  2  using  several  program  options. 

In  T-ink  l,  the  input  data  array  is  always  printed  out  in  matrix  form  as 
well  as  are  the  generated  arrays  characterizing  the  boundary  and  inter¬ 
face  points.  Upon  an  input  signal,  a  descriptive  printout  is  given  of 
the  allocated  data  array  and  the  remaining  generated  parameters. 

In  Link  2,  at  regularly  spaced  time  intervals  (the  spacing  may  change 
with  each  change  in  Fon(T)),  tke  following  computed  values  are  printed 


Pig.  8  Link  2  Procedure 
Plow  Diagram 


(a)  Surface  data  at  each  intersection  of  the  hot  gas  boundary  with  a 
radial  wash  .JLne,  including  mass  flowrate  of  gases  generated  within  the 
ablative  region,  new  positions  of  the  surface,  and  surface  recession 
rates* 

(b)  Temperatures  at  all  points,  regular  and  irregular,  which  lie  on  a 

horizontal  (vertical)  mesh  line  In  an  odd  (even)  time  step* 
x 

She  following  information  is  available  for  optional  printout  at  the 
same  times: 

(a)  Values  of  effective  heat  transfer  coefficients  at  both  ends  of 
each  horizontal  (vertical)  mesh  line  in  .an  odd  (even)  time  step* 

(b)  Temperatures  at  irregular  boundary  and  interface  points  lying  on 
vertical  (horizontal)  mesh  lines  in  odd  (even)  time  steps. 

(c)  Values  of  0  (x, y,  t)  and  R(x,y,*r)  at  every  regular  point  of  the 

y 

charring  ablator  which  has  reached  the  minimum  pyrolysis  temperature. 

Finally,  a  listing  is  printed  out  at  the  end  of  the  run  of  the  maximum 
temperature  achieved  at  each  regular  point  of  the  charring  ablator* 
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RESULTS 


la  order  to  provide  a  broad  base  for  evaluating  the  capabilities  of  the 
two-dimens  ional  computer  program,  2D-ABLATB,  developed  in  Phase  II  of 
this  study,  checkout  of  the  program  Included  comparisons  both  with 
theoretical  and  with  test  data,  as  follows: 

1*  Simulation  of  transient  one-dimensional  (radial)  conduction  in  a 
hollow  cylinder,  with  no  charring  or  recession.  The  exact  solution  for 
this  problem  is  known  and  is  presented  in  the  form  of  temperature 
response  curves  in  Bef  •  21,  the  "CYLEEAT"  handbook. 

2.  Comparison  with  test  data  obtained  from  an  engine  firing  at  Edwards 
APB.  Here  it  was  necessary  to  simulate  two-dimensional  charring  and 
conduction  in  a  two-material  thrust  chamber,  the  exposed  material  being 
a  carbon  cloth-phenolic  backed  up  by  a  stainless  steel  shell. 

Cases  1  and  2  above  comprised  only  a  portion  of  the  checkout  performed 
for  the  2D-ABLAXB  program,  which  for  purposes  of  debugging  Included 
treatment  of  an  assortment  of  wall  configurations,  material  arrangements 
and  properties,  and  external  heating  conditions.  Cases  1  and  2,  however, 
were  treated  in  depth  and,  being  representative  of  the  two  classes  of 
problems  handled,  are  discussed  here  in  detail. 


O 
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COMPARISON  WITH  AN  EXACT  SOLUTION  OF  A  CONDUCTION  PROBLEM 


O 


Despite  the  complexity  of  the  2D- ABLATE  program  and  the  many  physical 
and  chemical  mechanisms  it  treats,  its  basic  function,  nevertheless,  is 
the  prediction  of  transient  temperature  distributions  resulting  from  the 
conduction  of  heat  in  solids  •  Consequently,  before  the  program  could  be 
used  with  any  real  degree  of  confidence,  it  vas  necessary  to  determine 
its  capabilities  for  handling  pure  conduction  problems  •  Furthermore, 
because  an  implicit  method  of  solution  vas  employed  for  the  energy  equa¬ 
tion,  it  vas  especially  important  to  show  that  sufficiently  accurate  and 
stable  results  can  be  obtained  using  large  distance  increments  and  rela¬ 
tively  large  time  steps,  l.e.,  lor&e  (for  a  given  distance  increment) 
compared  to  the  step  size  required  for  an  explicit  method.  For  this 
purpose,  the  best  standard  for  comparison  is  an  exact  solution.  Although 
the  case  chosen  of  a  hollow  cylinder  is  only  one-dimensional,  the  spatial 
coordinate  is  radial  (the  direction  of  greatest  heat  flux  in  most  engine 
firings)  and  provides  the  closest  approximation  to  an  axisyanetric 
thrust  chamber  for  which  an  exact  solution  is  available. 


The  conduction  problem  whose  solution  is  given  in  the  CYLHSAT  handbook 


can  be  expressed  (in  our  notation),  'ab  follows: 
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r  &r  n 

y^o 

(69) 

T.T„ 

» 

T  ■  0 

(70) 

where  a,  K,  and  h  are  known  constants  and  where  r .  and  r  are  the 

1  o 

Inner  and  outer  radii,  respectively,  of  the  cylinder  (see  Pig.  9  for 
axial  and  radial  cross-sections) . 

In  order  to  assess  the  2D-ABJATE  conduction  capabilities  under  extreme 
conditions,  the  program  was  run  for  the  cases  of  both  high  and  low 
thermal  conductivity  materials,  such  as  an  ATJ  graphite  throat  insert 
and  a  carbon  cloth-phenolic  ablator,  respectively.  The  numerical  data 
used  as  input  for  the  two  cases  are  given  in  Table  3.  Only  a  single  list 
of  data  is  given  to  cover  both  cases  because  the  values  chosen  differ 
only  In  thermal  conductivity.  The  numerical  solution  was  obtained  in 
both  cases  using  a  distance  increment  of  .11  inches.  Multiple  runs 

were  made  in  each  case  for  several  different  constant  values  of  At, 
ranging  from  .125  ts  2  seconds  in  Case  1  (the  high  conductivity  case)  and 
from  .5  to  8  seconds  In  Case  2.  The  results  are  presented  (for  At  *  1 


Table  3«  Properties  and  Conditions  Used  for  the  Hollow  Cylinder  Runs 


ami  8  seconds  only  in  Case  1  and  for  At  »  2  and  3  seconds  in  Case  2)  in 
Figs.  10  to  13  in  the  forts  of  surface  temperature  histories  and  tempera¬ 
ture  profiles  at  fixed  tine  levels. 

The  power  of  using  a a  top  licit  method  is  strikingly  illu  Jrated  here, 
especially  in  the  high  conductivity  case.  Had  ve  used  the  explicit 
method,  ve  see  from  Inequality  (65)  that;  to  ensure  a  stable  solution, 
ve  vould  have  been  forced  to  use  a  value  of  At  less  than  or  equal  to 
.1J2  seconds  in  Case  1  and  .906  seconds  in  Case  2.  Tot,  as  shown  in 
Figs.  10  and  12,  ve  vere  able  to  obtain  a  stable  solution  (and  quite 
accurate  -  especially  at  the  surface)  using  a  tine  step  as  nnch  as  60 
tines  as  large*  Theoretically  ve  can  use  as  large  a  tine  step  as  ve 
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FIG.  12.  HOLLOW  CYL1MX3R 

of  Interior  at  64-  seconds  -  High  Thermal  Diffusivity  Material 
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2Siaraal  Response  of  Interior  at  60  seconds  -  Low  Theraal  Diffusivity  Material 


please  and  still  expect  a  stable  convergent*  solution  -  it  would  Just 
take  longer  to  home  in. 

In  Fig*  10,  for  both  At  ■  1  second  and  hi  =  8  seconds,  ve  see  that  the 
calculated  surface  temperatures  begin  below  the  exact  curve  and  approach 
it  asymptotically  with  time*  Simultaneously,  a  small  error  is  evident 
in  the  numerical  solution  in  the  form  of  a  slight  oscillation  which 
damps  with  time*  The  larger  one-sided  error  is  truncation  error  due  to 
discretization  of  the  continuous  Fq.  (67},  whereas  the  oscillatory  behav¬ 
ior  is  characteristic  of  the  numerical  procedure  itself  independent  of 
the  differential  equation.  In  particular,  the  damping  of  the  oscilla¬ 
tions  with  time  illustrates  the  stability  of  the  procedure*  The  computer 
run  made  with  hi  *  1  second  was  terminated  after  8  seconds  due  to  the 
expected  early  convergence  of  the  numerical  with  the  exact  solution.  For 
the  run  performed  with  so  large  a  time  step  as  8  seconds,  on  the  other 
hand,  it  was  considered  advisable  to  carry  out  the  computations  for 
almost  1J0  seconds  to  ensure  convergence*  As  can  be  seen  in  Fig*  10, 
this  proved  to  be  overcautious;  for,  despite  an  error  of  over  500°R  at 
16  seconds  (the  first  available  surface  calculation),  most  of  the  early 
error  has  been  eliminated  shortly  after  50  seconds  and  the  correspondence 

♦Although  the  criteria  given  above  for  size  of  the  time  step  were  derived 
in  terms  of  stability  of  the  numerical  procedure,  it  can  be  shown  that 
this  is  intimately  related  to  convergence*  Indeed,  In  the  case  of  a 
linear  problem  such  as  that  defined  by  Eqs*  (67)  -  (70),  a  proof  is  given 
in  Ref*  22  that  stability  and  convergence  are  equivalent  for  a  consistent 
numerical  scheme,  i.e*,  one  for  which  the  difference  equation  at  each 
point  approaches  the  differential  equation  in  the  limit  as  the  space  and 
time  Increments  go  to  sero* 


is  good  considering  the  size  of  the  time  and  space  increments. 

The  numerical  results  obtained  for  the  low  conductivity  material  are 
pictured  in  Fig.  11.  The  closer  approximation  to  the  exact  solution  is 
to  be  expected  because  an  8  second  time  step  here  is  not  as  large  rela¬ 
tive  to  that  required  for  the  explicit  method  as  it  vas  in  the  high  con¬ 
ductivity  case.  The  reason  for  the  reversed  directive  of  the  early  error 
in  the  2  second  time  step  run  is  not  evident* 

Figures  12  and  13  compare  temperature  profiles  computed  by  2D-ABLASQ2  for 
At  «  8  seconds  with  the  exact  solution  at  fixed  representative  inter¬ 
mediate  time  levels  selected  late  enough  to  avoid  early  gross  truncation 
error  hut  before  most  of  the  error  has  been  damped  out.  Again,  the 
errors,  although  not  excessive,  are  greater  in  the  case  of  the  high  con¬ 
ductivity  material.  In  particular,  the  error  can  be  seen  to  increase 
with  depth  into  the  material  in  the  high  conductivity  case,  whereas  this 
behavior  is  not  particularly  evident  for  the  lav  conductivity  material. 
This  should  not  be  too  unexpected  because,  with  increased  distance  from 
the  surface  exposed  to  heating,  the  stability  (and  thus  convergence) 
criterion  given  by  Inequality  (63)  becomes  less  applicable  and  is  re¬ 
placed  by  Inequality  (28),  which  covers  conduction  in  the  Interior  unaf¬ 
fected  by  external  heating.  On  the  basis  of  Inequality  (28),  the  upper 
bound  for  &r  required  ixSSthe  explicit  method  remains  virtually 


! 


unchanged  for  the  high  conductivity  material  but  increases  significantly 
for  the  low  conductivity  material,  taking  on  values  of  .151  seconds  and 
7*56  seconds,  respectively.  Thus  ve  see  that,  for  internal  conduction, 
a  time  step  of  8  seconds  for  the  high  conductivity  material  is  still 
more  than  50  tines  as  large  as  that  required  for  the  explicit  method, 
but  is  only  slightly  larger  for  the  lew  conductivity  material.  Thus  ve 
expect  much  more  accurate  results  in  the  latter  case.  Indeed,  for  a 
low  conductivity  material,  aside  frm  reduced  local  temporal  truncation 
error,  there  is  little  advantage  in  the  use  of  an  implicit  method  as 

long  as  the  spatial  temperature  gradients  are  not  too  steep.  When  large 

/ 

gradients  occur,  however,  small  distance  increments  are  required  for 
local  spatial  accuracy  and  an  implicit  method  vould  permit  correspond¬ 
ingly  larger  time  steps  when  desirable  than  vould  the  explicit. 

1 

In  summary,  Figs.  10  through  13  show  that  2D- ABLATE  can  provide  very 

accurate  results  for  a  purs  conduction  problem  using  relatively  large 

distance  and  time  increments.  Increased  accuracy  vould  undoubtedly 

result  from  exercising  the  program  option  of  varying  the  time  step 
# 

during  a  computer  run,  using  a  smll  step  at  early  times  and  larger 
steps  later. 


COMPARISON  WITH  ABLATIVE  ENGINE  TEST  MTA 


After  establishing  the  adequacy  of  the  cB-ABLAIE  progiaa  in  perforating 
its  basic  function,  i*e v,  in  the  treatment  of  heat  conduction,  the  next 
step  was  to  evaluate  its  usefulness  in  practice*  For  a  real  engine 
firing,  the  program  must  be  able  to  predict  transient  temperature  dis¬ 
tributions  with  a  reasonable  degree  of  accuracy  within  the  limits  of 
a}  the  aunlicability  of  the  mechanisms  simulated  in  the  program  n-!^ 
b)  the  accuracy  of  the  values  used  as  input  for  the  physical  and  chemi¬ 
cal  properties*  Furthermore,  in  order  to  justify  the  use  of  an  implicit 
scheme  for  the  numerical  solution,  the  time  steps  employed  in  the  pro¬ 
gram  must  be  larger  than  those  required  for  the  explicit  method*  To 
demonstrate  this  capability,  a  case  vas  run  on  2D-ABLATE  for  which  test 
data  were  available  (supplied  by  Edvards  Air  Force  Base)  •  Reasonably 
good  agreement  with  the  measured  temperature b  and  char  depths  was  ob¬ 
tained,  at  least  in  the  threat  region,  as  is  shown  below*  The  data 
supplied  were  for  n  thrust  chamber  firing  Fg/HgH^  blend  (66.7$ 

24jt  MMH,  and  9*5#  HgO,  by  weight)  at  a  chamber  pressure  of  116  psla  to 
develop  a  thrust  of  5900  pounds*  The  thrust  chamber  vails  consisted  of 
two  materials  including  a  phenolic-carbon  cloth  ablative,  exposed  to 
engine  firing,  backed  up  by  a  stainless  steel  shell  (see  Fig.  14) *  The 
external  heating  (and  cooling)  conditions  included  heating  of  the  inner 
surface  by  the  hot  combustion  gases,  radiative  cooling  at  the  outer 


Gurfaco  and  at  the  diverging  portion  of  the  inner  nozzle  surface,  con¬ 
vective  cooling  (faimulatcd  in  the  radiation  term)  at  the  outer  surface, 
end  insulation  elsewhere  as  indicated  in  Fig*  14*  Engine  firing  vas 
steady  for  the  first  60  seconds,  fallowed  bv  edglnt  shutdown  and  soakout 
for  the  remainder  of  the  run*  The  a^daJTheat  flux  in  the  combustion 
chanter  vail  vas  assumed  to  be^so^slight  that  an  effective  insulated 
left  boundary  vas  esployed  at  an  axial  position  18  Inches  from  the 
injector  for  the  computer  run,  as  shown  by  the  dotted  line  in  Fig.  14; 
i*e«,  results  to  the  left  of  the  assumed  boundary  were  assumed  to  be 
Identical  to  those  obtained  to  the  right  (in  the  combustion  chamber 
only)  •  Sable  4  is  a  list  of  the  other  pertinent  data  required  to  des¬ 
cribe  the  firing  and  used  as  input  to  the  computer  program*  .Included 
are  material  and  gas  properties,  both  physical  and  chemical,  as  fcell  as 
the  operating  conditions,  heating  parameters,  erosion  date,  and  the 
space  and  time  increments  used  to  obtain  a  finite  difference  solution* 

The  values  used  vere  based  on  information  supplied  by  EAFB  supplemented 
by  additional  Rocketdyna  data  and  experience*  The  density-specific  heat 
product  employed  for  the  ablative  material  accounts  for  both  the  transi¬ 
tion  from  virgin  material  to  char  and  the  increase  of  the  specific  heat 
of  the  char  with  tesperature*  The  thermal  conductivity  of  the  ablative 
reflects  similar  effects  in  addition  to  a  radiation  contribution  at 
high  temperatures*  Pyrolysis  of  the  ablative  resin  occurs  over  a 
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TfflJS  4.  FEOPSKTISS  ABB  COBBmDHB  USED  TO  SBCJIATE  REAL  BtOBB  mOB 


HA3SRIAL  PROPERTIES 


Phcnollc-Carbon  Cloth  Ablative  (and  qenerftted 


ysla  ga ub): 


PC  -  I.38XIO"2  3tu/in3-°R  ,  460°  sis  o60°a  , 

.42XlO*5T  ♦  .9768X10"2  ,  960°  <TS  i9600  , 

,5S10“35?  lt212XlO"2  >  1960°  SI  ^  5460°  , 


2.25XL0'2 

,  3460°  sc  T  *s  6960°  • 

K  -  •8X10"8T  +  .7J2X10-5 

Btu/in-aec-°R  ,  460°  *  T  <  96C°R  , 

l.XlO^T  +  .54xl0“5 

,  960°  sts  3160°  , 

,5X10“^T  +  2 .12X10"5 

,  3160°  sis  696O0  • 

H«  .5T  -  .023  Btu/lba 

,  460°  sis  696o°r  . 

*n.  -  »5 

F  a  ° 

py 

,  is  960°  , 

>5  [i.  SS-^21] 

,  960°  <T<  I9600  , 

•6 

Qpy  a  25O  Btu/lb  resin 

,  I9600  *5  T 

P*  -  .0542  lba/in5 


TABIS  4  (Coat'd.) 


© 


504  Stainless  Steel: 

PC  -  .031  Btu/in5-°R  ,  460°  ST  s  6960 °R 

K  -  1.X10"7T  +  1.54x10^  Btu/ln-s8c-°R  ,  460°  s  T  s  6960°B 

EROSION  CONSTANTS  AND  COMBUSTION  GAS  CONDITIONS 
(See  NonenaLature  in  Appendix  c) 

Gas  Plow  Constants: 

7  -  1.312 

-  18.574  V“We 

Q  Pc  »  116  psia 

(Cp)w  »  .452  Btu/lbm-°R 
**c  ■  3.2335  in. 

vaporisation  Constants: 

D  m  ,02  ln2/sec 
vap  • 

-  -Si06*105  SWU'.OK, 

»  4000  Btu/lbm 

Pvmll  Material  “  *°^5 


© 
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IAELS  4  (Cont'd.) 


Reaction  Consents  ( C ^  -*  ®4(g)  +  C2HE(g)^! 


'reactant  gas 


2  lb  /lb 

mr  i 


nr  mole 


1 reactant  gas  “  ^ 

»  .0089  in/sec 

w  ■  ^ 


reaction 


5000  Btu/lb 


T\  y 

reaction  »  .118  in  /sec 

(cpW=u  -  -3  BWV* 


(SESAHnSG  CONDITIOSS  AS)  HEATI5G  EABAMSEEBS  (axially  varying) 


Tc  "  53<>OR 

T  » °°R 
env 

B  =  .4 


•452  Btu/lbffl-°R 

i\py  ■»  1.  (for  pyrolysis  gases  ejected  at  surface) 
ob  »  3.3  X  10~15  Btu/in2-sec-(°R)^ 


TAELS  4  (Coat'd.) 


TABIS  4  (Coat’d.) 


FIHTE3  DIFFERENCE  INFORMATION 

Space  Increments* 

y  m  0,  at  axis  of  thrust  chamber 

by  m  .I  An,  axial  mesh  lines  at  y  a  ■  3,3  +  (j-1  )isyt  j  *  1,  • 

x  ®  ©,  ,1  in.  to  left  of  chaster  out-off  position 

radial  aesh  lines  at  x  a  xi+1  -  +  o*^  i  ■  I,  ‘**,14,  x^  ®  o35 


i  12345 

6  7 

8 

9 

10  11  12 

15  14 

4xi  (in.)  .75  .75  .75  1.25  1.25  1*25  1. 

.  *5 

*25 

.25  .5  .75  1*5  1*1 

Variable  Time- Steps  (Option  A): 

k  1  2  3  4  5  6 

7 

8 

9 

10  11 

12 

AT^sec)  .1  .2  .5  1*  2.  5. 

1. 

2. 

5. 

10.  20. 

50. 

Fon,k^  1*  1#  1#  1*  1#  1# 

Co 

0. 

0. 

ff- 

O 

• 

O 

0. 

t'  0.  .3  .5  1.  10.  30. 

60.” 

70. 

80 

.  100.  200 

.  400. 

where  1  is  initial  time  for  use  of  4?^ 


teeperaturo  ran go  extending  froa  500  to  1500°F  vlth  a a  average  teat  o t 
pyrolysis  of  250  Btu/lb  4 In.  The  pyrolysis  gases  generated  were 

assused  in  this  ease  to  undergo  no  further  high  temperature  cracking  or 
gas-char  reactions  during  passage  through  the  porous  char*  The  gas- 
side  heat  transfer  coefficient  distribution  used  in  the  prograa  Is  based 
on  the  Barts  equation*  Thermodynamic  and  kinetic  constants  are  given  in 
Table  4 -to  account  for  possible  vaporisation  of  the  char  and  reaction  of 
the  char  vlth  hydrogen  in  the  coob^siion  gas*  As  noted,  *.  variable 
axial  distance  Increment,  dr,  was  used  with  closest  spacing  in  the 
throat*  A  variable  time  step,  At,  vas  specified  to  obtain  sufficient 
accuracy  with  a  minimum  of  c  cep  uter  time* 

The  pertinent  results  of  the  computer  run  of  the  2D-ABLATE  program,  for 
the  data  listed  in  Table  4,  are  presented  in  Figs*  15  to  19  in  the  form 
of  temperature  and  radial  char  depth  histories,  temperature  profiles, 
and  contours  of  the  char  front  at  fixed  time  levels*  Although  the  com¬ 
puter  run  included  erosion  constants  as  input  data,  no  surface  recession 
vas  predicted*  This  vas  in  agreement  vlth  the  results  of  the  engine 
test  firing* 

Figure  15  compares  computed  char  depth  histories  at  the  thre&t,  combus¬ 
tion  chamber,  and  exit  (at  axial  positions  4*87,  U .87,  and  .67  inches, 
respectively,  from  the  exit  plane)  with  tvo  experimental  threat  values 
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OS.  16.  SURFACE  TEMPERATURE  RESPOHSE 
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obtained  frca  ceasured  tost  data*  For  all  plotted  points  except  the 

neaoured  throat  value  at  700  seconds  (a  depth  of  1*2  inches),  char 

depth  for  a  given  axial  position  at  any  time  was  defined  to  be  the  radial 

distance  froa  the  surface  of  the  point  whose  aaxinua  temperature  attained 

up  until  that  tine  vas  1460°R,  the  midpoint  of  the  charring  temperature 

range  (the  60  second  experimental  value  vas  obtained  by  interpolation 

from  throat  temperatures  measured  at  several  radial  positions)  •  The 

measured  value  at  700  seconds,  on  the  other  hand,  vas  not  based  an  the 

location  of  a  temperature  level  but  instead  on  &  visual  Judgment  vhich 

placed  it  in  the  initial  portion  (considerably  below  1460°R)  of  the 

charring  temperature  range*  The  extrapolated  experimental  char  depth  at 

the  throat  (baaed  on  the  lh60^R  definition)  at  700  seconds  thus  should 

appear  much  closer  to  the  computed  value*  For  the  throat,  chamber,  and 

exit  values  of  T  and  £ _ employed  in  the  program,  the  relative 

av  conv 

positions  of  the  corresponding  computed  char  depth  histories  are  consis¬ 
tent  with  expected  results* 

In  Fig*  16,  surface  temperature  histories  are  presented  for  the  same 
three  axial  positions  (throat,  chamber,  and  exit)*  Although  no  measure¬ 
ments  were  made  of  surface  temperature  during  the  real  engine  firing, 
the  curves  shown  in  <ig*  l6  provide  further  grounds  for  evaluation  of 
the  run  based  on  l)  their  positions  relative  to  each  other,  2)  how 
closely  the  magnitudes  attained  correspond  to  expected  results,  and 


3)  the  compatibility  of  the  shapes  assumed  by  the  curves.  The  results 
sees  reasonable  on  the  basis  of  all  three  criteria  (in  assessing  the 
compatibility  of  the  shapes  assumed,  it  aust  be  remembered  that,  due  to 
axial  ordering,  the  ch&eber  and  exit  curves  should  not  be  compared 
directly  with  each  other,  but  only  through  comparisons  of  each  vith  the 
throat  curve). 

Figure  17  Is  a  comparison  vith  experimental  data  of  a  computed  tempera* 
ture  profile  at' the  throat  after  60  seconds  of  firing  (comparable  data 
at  the  chamber  and  exit  vere  not  available  in  depth) .  In  general,  the 
agreement  is  reasonably  good.  Ho  comparison  is  possible  until  a  depth 
of  *3  inches,  vhere  the  computed  value  is  50°R  less  than  the  measured. 
Proceeding  into  the  material,  the  computed  temperature  does  not  drop  off. 
as  quickly  and  is  above  the  next  measured  value  (at  a  depth  of  .73 
inches)  by  over  150°R.  Beyond  that  point,  there  is  not  enough  thermal 
penetration  for  a  meaningful  comparison. 

In  Figr  IS,  a  comparison  is  given  of  the  computed  temperature  history 
vith  the  measured  at  a  depth  of  .3  inches.  Again  the  correspondence  is 
reasonably  good.  Use  slopes  are  close  and  the  values  differ  by  at  most 
250%.  The  computed  curve  rises  sooner  than  the  measured,  the  differ¬ 
ence  between  the  curves  reaching  a  maximum  after  shoot  30  seconds  of 
firing.  At  this  tlsw,  the  slope  of  the  measured  curve  has  just  about 


c&u$it  up  with  and,  thenceforth,  exceeds  that  of  t)s  ccegputed  carve*  The 
curves  cross  at  about  55  seconds  and  the  computed  curve  has  dropped  below 
the  measured  by  about  60°R  at  60  seconds,  the  end  of  firing* 

Similar  comparisons  are  not  pictured  for  the  chamber  and  exit  tesperature 
histories  (no  similar  test  data  are  available  at  the  exit)*  Ths  cor¬ 
respondence  in  the  chamber  vas  not  good,  there  being  almost  a  JO  second 
delay  in  the  beginning  of  the  rise  of  ths  measured  temperature  values 
after  that  of  the  computed  values*  Xo  physical  explanation  has  been 
found  for  this  lag  unless  some  other  heat  transfer  mechanism,  such  as 
film  cooling  from  the  injector,  vas  present  in  the  combustion  chaaber  to 
cause  the  delay*  The  computed  chamber  values,  on  the  basis  of  past 
experience,  to  be  more  compatible  with  both  the  computed  and 
measured  throat  .values .  After  the  Initial  lag,  the  slope  of  the  curve 
of  the  measured  values  is  relatively  steep  end,  by  the  end  of  the  60 
ascend  firing  cycle,  the  ssasured  value  is  within  50°B  of  the  computed • 

Finally,  in  Fig*  19,  contours  are  pictured  of  the  char  front  (based  on  a 
146C°R  isotherm)  obtained  from  the  computer  run  for  two  fixed  time  levels 
at  *  »  60  seconds,  the  end  of  firing,  and  at  t  ■  700  seconds,  after 
640  seconds  of  soakout*  It  is  of  interest  to  note  that,  although  the 
radial  char  depth  at  60  seconds  is  nearly  uniform  for  all  axial  positions 
(see  Fig*  1$  as  veil),  the  additional  radial  penetration  of  the  char 


front  during  sookout  is  considerably  greater  at  the  throat,  then  else¬ 
where,  the  least  additional  penetration  occurring  at  the  exit  position* 
This  is  mainly  due  to  the  greater  heat  content  in  the  throat  volume 
below  the  char  front  at  60  seconds  than  in  the  exit  or  cfcassber*  In 
addition,  the  particularly  snail  change  at  the  exit  is  probably  due 
in  part  to  axial  conduction  losses  from  the  exit  to  the  throat  plane  and 
possibly  to  greater  proximity  of  the  exit  portion  (and  the  chamber  por¬ 
tion  as  well)  of  the  char  front  to  the  radiating  sink  at  the  outer  sur¬ 
face  of  the  stainless  steel  shell*  in  the  first-  60  seconds,  on  the 
other  hand,  these  effects  are  dominated  by  the  strong,  almost  uniform 
heat  flux  along  the  exposed  Inner  'vurface*  Although  no  contours  of 
measured  data  are  available  for  comparison,  the  computed  contours  agree 
with  the  60  eecond  experimental  point  pictured  in  Fig*  1?,  they  assume 
shapes  which  are  in  agreement  with  the  relative  locations  of  the  adja¬ 
cent  thermal  sources  and  sinks,  and  they  are  compatible  with  each  other* 

To  suaaarize,  the  results  of  using  the  21VABULT3  program  to  simulate  a 
real  engine  firing  corresponded  reasonably  well  to  measured  test  data, 
and,  where  measured  data  were  not  available,  the  computed  results  v«re 
compatible  with  expected  results  and  with  each  other.  They  were  obtained 
for  a  TOO  second  run  using  time  steps  which  varied  from  .1  to  30  seconds 
in  size  (see  Table  4),  79  steps  in  all*  The  machine  time  required  for 
the  run  can  he  broken  down  as  follows;  30  seconds  to  load,  20  seconds 


■i^i^'na.  -  «<jri<syga!ta6Siaa8S> 


for  the  Link  1  preliaicary  execution,  and  120  seconds  for  the  tisss-step 
calculation  (at  seconds  of  ss&chlne  tine  per  tine  step),  for  a  total 
of  170  seconds  of  machine  tiae.  On  the  basis  of  Inequality  (28)  (implied, 
to  the  backup  material,  dv»  to  its  high  tfeara&l  diffusivity)  the  expli¬ 
cit  method  would  have  required  a  value  of  At  less  than  or  e>ssal  to 
,5  seconds.  Even  assuming  as  little  m  §  second.  0?  mehine  time  per 
time  step,  this  would  scan  a  jainism  of  TOO  seconds  for  execution  alone 
wiiSo  ^  explicit  setheds  Furthexsoxe,  t-****  nuebsr  of  tlsss  stsps  sod 
thus  the  raehine  tine  us©&  by  2D«fi8yws  for  this  caso  could  probably 
have  been  reduced  significantly  without  appreciable  changes  in  the 
results.  This  was  substantiated  by  making  another  7§0  second  run  Men~ 
tical  to  the  first,  but  for  which  each  the  step  was  cut  in  half*  The 
temperatures  obtained  at  ccsQparable  tises  and  positions  were  uniformly 

within  ljt  of  each  other  from  run  to  run  (except  in  the  first  few  steps 

\ 

of  both  the  firing  and  cooling  periods  when  they  differed  by  at  nest  yf») 
and  generally  differed  by  considerably  less  than  1 £•  This  is  veil 
within  the  usual  convergence  standards  and  indicates  that  larger  steps 
than  those  used  in  the  fl-'st  run  could  have  been  safely  employed. 
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CONCLUSIONS 


In  Phase  2  of  the  subject  program,  It  has  been  demonstrated  by  means  of 
a  digital  computer  code  that  the  nonlinear  equations  of  a  comprehensive 
two-dimensional  transient  ablation  (including  both  surface  erosion  and 
internal  charring)  and  heat  conduction  model  for  the  multimaterial  walls 
of  a  rocket  engine  thrust  chamber  can  be  solved  using  numerically  stable; 
implicit  finite  difference  procedures.  The  analysis  was  performed  in 
cylindrical  coordinates  (axial  and  radial)  assuming  axisyisaetrie  mils 
with  temperature  dependent  physical  and  chemical  properties  and  arbi¬ 
trarily  curved  material  boundaries  and  interfaces.  Chemical  reactions 
within  the  charring  material  vere  treated  in  depth  and  include  three  gas 
generation  reactions  plus  cracking  of  the  generated  gases.  Surface  ero¬ 
sion  effects  handled  at  the  hot  combustion  gas  boundary  vere  malting, 
vaporization,  and  chemical  reactions.  Blocking  of  the  combustion  gas 
due  to  injection  of  gases  generated  within  the  vail  materials  and  at  the 
surface  was  also  accounted  for.  Other  heat  flux  mechanisms  included  in 
the  model  and  handled  by  the  program  included;  engine  firing,  environ¬ 
mental  heating,  radiation  from  the  vails  to  the  environment,  insulation 
at  the  boundaries,  and  convective  cooling  within  the  charring  ablator. 
Engine  firing  say  be  steady  or  intermittent  with  sookout,  the  latter 
consisting  of  either  multiple-start  operation  or  pulsing.  All  heat 


flaxes  et  the  boundaries  and  Interfaces  were  expressed  In  terms  of 
normally  directed  temperature  gradients  . 

A  tine  step  procedure  was  employed  to  uncouple  and  solve  the  continuous 
equations,  finite  differences  being  used  for  the  discretization.  Of 
particular  interest  was  a  generalization  of  the  unconditionally  stable 
alternating  method  of  Peaceman  and  Rachford  used  for  the  energy  equation 
in  conjunction  with  special  second  order  accurate  techniques  developed  * 
to  simulate  the  normal  gradient  conditions  at  curved  surfaces «  The 
numerical  procedures  developed  were  incorporated  in  the  Fortran  IV  com¬ 
puter  program,  2D-ABLATE,  used  for  the  solution. 

The  program  capabilities  were  evaluated  by  comparison  with,  l)  an  exact 
solution  of  the  transient  one-dlmenslonal  pure  conduction  problem  for  a 
hollow  cylinder  heated  at  the  inner  surface  and  insulated  at  the  cuter > 
and  2)  test  data  from  an  engine  firing  with  two<  •dimensional  transient 
charring  and  conduction  taking  place  In  the  two-material  thrust  chamber. 
CcEpsrisona  were  mad'?  at  various  axial  positions  with  respect  to  char 
depth  (in  tLa  second  case  only),  temperature  profiles,  and  both  surface 
and  interior  t®ss>erature  histories.  Excellent  agreement  was  obtained 
with  results  of  the  exact  solution  and  reasonably  good  agreement*  with 
the  measured  test  data  obtained  from  the  engine  firing.  Moreover,  the 
program  exhibits  a  high  degree  of  self-consistency  and  compatibility  of 
results  obtained. 
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RECOMMENDATIONS  FOR  FUTURE  EFFORTS 

The  Phase  1  effort  of  the  subject  program  has  led  tc  the  development  of 
numerical  procedures  and  a  computer  code  for  solving  comprehensive  two- 
dimensional  transient  ablation  and1  heat  conduction  problems.  There  are, 
however,  several  extensions  to  the  program  which  should  be  made  in 
addition  to  both  practical  and  theoretical  investigations  to  de¬ 
limit  the  program  capabilities  and  provide  analytical  tools  for  future 
rocket  nozzle  designs.  A  summary  of  recommendations  for  future  effort 
follows. 

1.  Many  materials  of  current  interest  are  strongly  anisotropic.  The 
existing  computer  program  should  be  extended  to  handle  materials  with 
differing  temperature  dependent  conductivities  in  the  x  and  y  direc¬ 
tions. 

20  Provision  should  be  made  in  the  program  to  permit  wzri-  than  one  of 
the  wall  materials  to  be  charring  ablators. 

3*  It  has  been  well-established  that,  for  many  cases  of  engine  firing, 
radiative  exchange  between  opposing  surfaces  of  the  hot  gas  boundary  has 
a  significant  effect  on  the  transient  temperature  distribution.  T^e 
program  should  be  extended  to  account  for  the  effect  of  this  mechanism. 

4.  A  comprehensive  study  should  be  made  of  several  of  the  parameters 
In  the  analysis  to  ascertain  their  significance  in  the  overall  heat 
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transfer  process.  The  goal  of  such  a  study  would  be  the  generation  of 
graphs  and  tables  which  summarize  the  variation  of  thermal  penetration, 
char  depth,  and  surface  erosion  in  rocket  engine  walls  with  the  essential 
parameters  governing  the  material  properties  and  arrangement  as  well  as 
the  thermal  environment.  Such  results  would  provide  basic  information 
required  for  the  design  of  thrust  chambers  and  would  serve  t'.  establish 
detailed  guidelines  for  the  subsequent  use  of  the  computer  program. 

5o  Investigation  should  be  made  to  develop  alternative  finite-difference 
procedures  which  provide  greater  accuracy  by  reducing  both  tern;  -ral  and 
spatial  local  truncation  errors.  The  principal  motivation  of  an  improved 
scheme  would  be  its  ability  to  reduce  the  machine  time  required  to 
achieve  a  given  standard  of  accuracy  by  increasing  the  size  of  the 
computational  increments  required.  A  three-level  (in  time)  alternating 
^wstion  method,  for  example,  which  has  recently  been  developed  by  the 
pyv'  snt  investigators,  might  be  employed.  This  procedure  is  second 
order  in  time,  requiring  no  linearization  for  its  application  to  the 
nonlinear  energy  equation.  Unconditional  stability  has  already  been 
established  for  this  method  when  used  to  discretize  the  linear  energy 
equation. 

60  A  theoretical  investigation  Bhovu •'  be  made  of  the  finite  difference 
procedures  utilized  in  tk'  as.. .  of  5  ts  extensions  ur  modi¬ 

fications  (in  accordance  with  v ,• <•  1..  2,  2,  and  5  above)  with 


main  emphasis  given  to  a  study  of  their  numerical  stability  and  con¬ 
vergence  properties.  The  value  of  such  an  investigation  would  be  as  a 
tool  to  r  xminate  expensive  experimentation  vith  the  program  to  establish 
detailed  quantitative  stability  and  convergence  criteria  for  each  change 
in  materials  and  thermal  environment.  In  the  event,  however,  that  a 
theoretical  investigation  is  not  made,  then  an  empirical  determination 
of  stability  and  convergence  criteria  should  be  conducted. 
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AEEENDK  A 


DEVELOPMENT  OF  ENERGI  AND  CONTINUITr  EQUATIONS 


The  following  form  of  the  maos  continuity  equation  has  heen  employed 
In  the  char  region  of  the  charring  ablative  material* 


(A-l) 


Using  the  theg&adynamlc  methods  .given  In  Appenduc  B  to  egress  the 
density  p  .of  the  charring  material,  ve  can  rewrite  Eg.  (A-l)  as 
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By  assuming  the  tin®  rate  of  change  of  gas  density  to  he  .dominated  hy 

f 

the  other  terms,  ve  obtain  the  form  of  the  continuity  equation  given 
earlier  {Eg. 

}fykl££  use  of  El*  (A~2%  -vs  can  further  derive  the  form  of  the  energy 
.equation  gftnsn  previously  hy  Eg.*  {!)*  We  begin  by  assuming  a  form 
for  the  (energy  agnation,  as  follows  {again  using  the  thermodynamic 
approaches 
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2!faen,  :multij>Iyin3  gg.*  (A-2)  by  E  and  adding  Sg,.  (A-3)  with  enthalpy 
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derivatives  expanded,  we  obtain  after  cpgfoinlng  terns  ajjprqpri^t^ty 
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Jly  rioting  that  0_  is  dominated  by  o,  we  can  ai®preso  the  second. 
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term  on'  the  left  aide  of  £&*  <V(A«$.*  3En  addition,  the  first  term  oan  Hie 
rewritten  as  foUcwss 
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where  .C  is  the  specific  heat  sat  constant  pressure  of  the  solid  s^te- 
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rial  in  the  ohar  region,  hie  xzy  similarly  rewrite  :US  3j5T:®f  * 
•Finally  we  aaay  take  to  he  o  constant  3^  by  nasnming  that 

the  ccEpomtion  of  the  ga3  generated  in  oach  reaction  remains  one  hanged 
os  it  is  being  f onaedi  i,s,,  for  -each  % 
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This  1.  s  by  integration  to  the  following  expression 
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Then  v  h|r  write  Eq,.  (A*4)  with  the  above  modifications 
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*.hich  i  Exactly  Eq.,  ( 1 ) .  (The  subscript  g  has  been  dropped  for 
venienc  from  the  symbol  for  gas  enthalpy.) 
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APPENDIX  B 


DERIVATION  OP  EQUATION  FOR  DENSITY  OP  CHARRING  LAMINATE 


The  derivation  of  the  egression  given  previously  by  Eq.  (  4 )  for  the 
density  of  the  solid  portion  of  the  charring  material  is  based  on  a 
thermodynamic  simulation  of  all  internally  occurring  chemical  reactions, 
and  in  particular  on  the  assumption  that  all  chemically  reacting  systems 
within  the  material  are  in  thermodynamic  equilibrium.  Required  for  the 
development  are  the  following  Known  functions  of  temperature: 

Fpy(TJ  -  fraction  of  mass  converted  from  resin  to  vapor  by  the 
pyrolysis  reaction  (then  l-F^  would  be  the  fraction 
converted  to  char,  i.e.,  the  fraction  of  solid  mass 
remaining  after  pyrolysis) 

F  (t)  -  fraction  of  mass  (of  char  and  reinforcement  combined) 

SC 

remaining  as  the  solid  product  of  the  char-reinforcement 
reaction 

Fdec(T)  -  fraction  of  mass  converted  from  the  solid  product  (of  the 

char-reinforcement  reaction)  to  gas  by  further  decomposition. 

Other  symbols  used  are  defined  in  the  Nomenclature  Section  of  this 
report. 

If  all  these  chemical  conversions  are  taken  into  account,  the  following 
expression  for  the  density  of  the  solid  material  for  any  degree  of  gas 


generation  results  from  a  mass  balance 


Since  each  F^  can  be  treated  as  a  ratio  of  densities  as  veil  as  a 
ratio  of  mass  increments,  ve  can  Immediately  vrite  by  definition 


O  S3  D*  p  S3  P*  P  P 

gl  res  py  v  res  py 
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The  density  of  the  intermediately  formed  gas  is  obtained  in  somewhat 
different  fashion  because  ve  ultimately  express  it  in  terms  of  the 
density  of  the  solid  product  instead  of  the  reactants.  Thus  we  first 
write 


P  .  3  P  +  P  -  O  . 

Hc  T  Hr  Hs 


(B-4) 


But  P  =  T _ P  and  P  =  F _ p  .  Therefore  we  obtain 

c  c,  s  s  r  r,  s  s 

pg8  ■  (ro,  s+rr,  *  «rc  s  +rr,  psc  • 

Finally,  substituting  Eqs.  (B-2),  (B-3),  and  (B-5)  in  (B-l)  yields 


P  -  P*(l-F  F_)  -  (r„  +T.  -l)p*  F  -  P*  F  F.  (B-6) 

y  res  py^  c,  s  r,  s  s  sc  s  sc  dec 


from  which  we  obtain  Eq.  (  4  )  by  combining  terms  • 


APPENDIX  C 


SURFACE  RECESSION  CALCULATIONS 


As  the  surface  temperature  of  the  thrust  chamber  wall  at  the  hot  gas 
boundary  increases,  several  modes  of  surface  removal  become  possible* 
They  include  melting,  vaporization  (or  decomposition),  and  chemical 
reactions  between  the  wall  material  and  certain  combustion  gas  compo¬ 
nents.  These  modes  of  surface  recession  were  considered  in  order  t.o  be 
able  to  handle  two  cases:  (l)  surface  removal  occurs  only  by  melting 
and/or  vaporization,  typical  of  glass-reinforced  ablative  materials; 

(2)  surface  removal  occurs  only  by  vaporization  and/or  chemical  reaction 
with  the  combustion  gas,  typical  of  carbon  or  graphite  cloth-reinforced 
ablative  materials  or  solid  graphite. 


For  both  cases,  the  energy  balance  at  the  hot  gas  boundary  can  be 
written  in  the  following  form,  equivalent  to  Eg..  (12)  (using 
a  -pv^  for  each  mode  i) : 
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Each  of  the  chemical  reaction  modes  is  governed  by  a  "kinetic  rate" 
equation  that  defines  the  temperature  dependence.  A  solution  for  the 
surface  recession  can  therefore  be  obtained  by  a  simultaneous  solution 
of  the  energy  and  kinetic  rate  equations.  A  one-dimensional  (radial 
coordinate  y),  constant-property  form  of  the  energy  equation  can  be 
written  in  the  form 


„  d2T  K  dT  _  v  dT  dT 
+  y  3y  **  Gg  -cp'g  oy  +  PC 


Integrating  Eq.  (C-2)  from  the  gas  surface  y  to  a  position  into  the 

8 

wall  material,  given  by  y  yields  (with  Eq.  (C-l)) 
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and  —  is  a  cons  ■feint  in  the  integration,  since  y  changes  only 
¥ 

slightly  in  the  Interval  (y  ,  y +s, )  •  To  evaluate  the  integral  of 

o  S  X 

Eq.  (C-3)  it  is  assumed  that  the  teaperature  profile  in  the  wall 


i  'J*  V,  S'  "t",  y  V  ^  y/-1 

v  s  5 

,/i  "  ■ 

?_ ... 


■  :,5Cc^-  r  -  5 .. 

-  -  f.  ■  ■  ~'S 


~'~V’*SS7*  *** 


*  *  y-  ^  -  ' 

'' '  v? 


material  can  be  approximated  by  the  quadratic  form, 
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T  a  T -a(y-y  )+b(y-y  )  ,  with  coefficients  determined  by  the  tempera- 
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tures  at  the  surface  and  the  first  and  second  mesh  points  on  each 
vertical  line.  Leibnitz* s  rule  gives 
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where  Svi  is  defined  by  Eq.  (C-^O)  or  (c-3l)  below.  Then  the  inte¬ 
gral  on  the  right  side  of  Eq.  (C-5)  can  be  obtained  as  follows: 

as,  be, 
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A  forward  time  difference  (which  we  abbreviate  here  using  the  operator 
8t)  is  vised  to  evaluate  the  time  derivative  giving 

ya+Si 
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Iteration  is  then  performed  on  the  value  of  surface  temperature  until 
Eq.  (09)  gives  the  same  recession  rate  Ev^  as  the  kinetic  rate 
equation  (C-27)  discussed  below. 


MELTING  ABLATION 

An  estimate  for  the  removal  of  material  by  melting  can  be  obtained  In 
simple  form  only  at  the  axial  position  for  which  the  surface  heating  is 
a  maximum  or  minimum  in  whose  neighborhood  the  surface  temperature  of 
the  melt  can  be  considered  constant. 


The  differential  equation  governing  the  flow  of  the  melt  layer  can  be 
written  in  the  following  form;  neglecting  the  effect  of  pressure 
gradients  (Ref.  17): 


where  6m  is  obtained  as  the  ratio  of  \(Es"Tm)  and  heff (Taw~Tg) . 
When  the  tangential  derivatives  of  heating  rate  are  zero,  £q.  (c-10) 


can  be  written  in  the  simpler  form, . 
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Reynolds  *  Analogy  can  be  used  to  relate  shear  stress  to  the  local  heat¬ 
ing  conditions: 
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In  the  throat  region  of  a  nozzle,  Sauer's  Analysis  (Ref*  18)  can  be 
used  to  estimate  the  velocity  gradient,  yielding 


3v"“ 


2  i1/2 

(r+l)rtrJ 


(C-14) 


In  addition,  the  melting  rate  equation  for  the  throat  can  be  written 
as 


v  -v 
v  vap 


con 


tli i.  r  «  t* 

^  e 


(c-15) 


For  purposes  of  the  present  analysis,  ve  assume  that  £g>  (c-ll)  holds 
for  the  entire  chamber  and  nozzle  regions  as  yell*  Finally,  the  free 
stream  velocity  gradient  is  estimated  on  the  basis  of  one-dimensional 
flew  to  yield 

^  ^ 2  ^  rt(l-.M*2} 
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For  use  in  £03.  (C-10)  -  (C-l6)  in  addition  to  the  equations  appearing 
below,  the  following  parameters  are  required: 


Pr*  JHL 

97-5 

(C-17) 

u  «  5  X  10"9  T?.*6 

5  w  oft 

(C=l8) 

(C-l8a) 
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(c)u 

“con/1 

(C-20) 
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(0-22) 
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MASS  TRAMSEER  EiFECTS 

In  computing  the  amount  of  surface  recession  for  reacting  or  valorizing 
materials,  the  diffusion  of  the  reactants  and/or  products  of  the  reac¬ 
tion  or  vaporization  through  the  boundary  layer  must  be  considered.* 

For  this  analysis,  a  frozen  composition  boundary  layer  is  assumed  where 
the  species  In  Question  diffuses  through  the  boundary  layer  without 
reaction.  In  this  case,  the  steady  state  maos  flux  of  the  species  Is 


governed  by  the  f oHovlng  differential  equation  (see  Bef  •  19)>  «b®»e 
Hie  subscript  1  is  emitted  for  convenience; 


p«d9 


-  p. 


Cc~23) 


Bor  a  constant  mass  diffusion  coefficient.  Eg.#  (C-23)  has  the  solution 

% 

*-*„+  <VV 

l-eD 

\  ' 

As  the  bloving  velocity  V  approaches  zero,  ve  obtain  the  following 
familiar  expression,  through  differentiation  of  Eg#  (C-24)  and  use  of 
lfHopitalfs  rule: 


V*- 

3E*1  ~T 

i 

i  v 

Shis  combined  with  a  particular  integral  of  Eg#  (C-23)  yields 


(C-23) 


pv  »  -p  y  *  - 
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peDW 

A 


(C-26) 


provided  the  solution  is  restricted  to  ¥«3 )/d  in  the  bouafiaay  layer# 


When  varorization  or  decomposition  occurs  as  the  mode  of  surface 
removal,  the  diffusing  canponcnt  is  considered  to  be  the  component 
vaporized  or  decomposed  at  the  surface#  Equation  is  used  to 

cccpute  the  recession  rate  because  the  mole  fraction  at  the  wall  is 
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M&xn,  as  discussed  la  the  asset  section.  Wien  &  chemical  reaction 
occurs  between  -the  eahaust  gases  and  the  wall  material,  the  diffusing 
component  Is  the  reactant  in  the  combustion  gas.  In  this  case  the  mole 
fraction  of  reactant  at  the  wan  must  he  computed. 


Hie  recession  rate  of  the  wall  Is  given  by 


-H/Jffl 


(C-2?) 


In  this  case  the  wall  concentration  Is  computed  from  the  following 
relation  obtained  from  £qs.  (C-26)  and  (C-2?) 


31  c  - 


therefore,  from  £qs*  (C-25)  and  (028)  we  obtain 


(C-28) 


^  t1  + — I  * 

© 

Equation  (C-28)  automatically  accounts  fear  diffusion  limited  reactions 
and  gives  a  smooth  transition  between  Mastic  controlled  and  diffusion 
controlled  react. ions. 


3h  Its  current  form,  the  computer  program  written  to  perform  -the 
surface  recession  calculation  can  ncccsscdata  a  rradmum  of 


three  Chemical  reactlo.us  occurring  IrdepeMently  at  the  material 
surface*  itotal  <eroalon  Is  (computed  as  thesmn  of  the  JrdlviduSl  ieac- 
tloa3  3>lu3  'Vaporization,  or  malting  plus  -valorization? 

Sv^  •*  Cohealcst  reactions)  ({C«50) 

or 

=■*  v^  {melting')  - 

^AK2B22»3!3m  CR  HECCKPDSHEICn 

Hhen  the  material  Sees  tot  melt,  anl  the  ethauct  gases  are  tou,-creact±ve 
to  the  parent  material,  snrtaee  remoral  is  ossnmed  to  occur  only  ly 
v^orisatlon  {in  the  ease  of  grag&dte)  or  thermal  decagatDgifcinn  {in  the 
oase  of  ^teflon,  silicon  carbide,  etc-),, 

‘Ihe  generation  of  vapor  is  handled  ms  a.  -^ero-crder  reaction  Iby  (com¬ 
puting  on  cg^Ilihriua  vapor  -pressure#  -dance  data  are  then  need  to 
(determine  a  vapor  pressure  sanation  <ef  the  tom 

-  :k  « 

vsp  yap 

Hhe  latent  isat  of  vaporization  (or  deco2?C3ltiaiO  is  then  ooual  to  33 
and  the  rate  of  surface  removal  (estimated  hy  diffusion  of  the  products 
of  decomposition  through  the  -boundary  dryer  113103  Xg-  {C-2S)# 
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1  Letters 

l  ' 

Specific  heat  of  vail  material 
'1  Specific  heat  at  constant  pressure 

*  f 

. «-  Coefficient  of  mass  diffusion 

*  ? 
Activation  energy  of  reaction 

Mass  flow  rate  per  unit  area 

Enthalpy  change  during  reaction 

Thermal  conductivity 

Mach  number 

I 

|  Molecular  weight 
|  .Pressure 

'  Prandtl  number 

I 

|  -  Universal  gas  constant 
*  Temperature 

i 

•  Adiabatic  wall  temperature 
Temperature  at  first  interior  mesh  point 
Temperature  at  second  interior  mesh  point 

i  Gar  velocity 

Blowing  velocity  of  ejected  gas  (®  -Pv/p) 
!  e 
|  Mole  fraction  of  gas 

* 

}; 
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Lower  Case  Letters 


cf 

f(x,T) 

bCORV 

heff 

\b& 

K 


V 

X 


y 


Sonic  Velocity- 
Coefficient  of  friction 
Radial  position  of  hot  gas  boundary 
Gravitational  constant 
Convective  heat  transfer  coefficient 
Effective  heat  transfer  coefficient 
Constant  factor  in  melt  viscosity  equation 
Constant  factor  la  vapor  pressure  equation 
Constant  factor  in  rate  equation 
normal  coordinate  to  surface 
Radius  of  curvature  of  the  throat 
Throat  radius  (value  of  y  at  -fee  throat) 

Normal  recession  rate  of  surface  (»  ^~  /  +  (df/dx)^) 

Axial  coordinate 
Radial  coordinate 


4 


Greek  tetters 


3  Angle  of  vail  surface 

y  Specific  heat  ratio  of  gases 

S„  Kelt  layer  thickness 

A  Effective  thickness  of  boundary  layer 

v  Tangential  coordinate  to  surface 


157 

/ 


/ 


{4  Viscosity 

P  Density 

f  Tijss 

Shear  stress 


a 

a 

g 

i 

m 

8 

V 

82 

vap 

v 

« 


Chamber  coalitions 

Conditions  '  at  edge  of  boundary  layer 

Conditions  in  gas  generated  vithih  chfcrrisig;  ablator 

Index  ranging  over  modes  of  vail  erosion. 

Melt  properties  or  soaditiona 

Gas-liquid  interface  (or  gas-vall  interface  if  no  liquid 
is  present)— often  referred,  to  os  '’surface”  in  te:ct 

Distance  from  surface  to  first  interior  mesh  point 

Distance  from  first  to  second  interior  mesh  point  j 

^spor  properties  or  conditions  f 

Wall  surface 

Free  stream  conditions 
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APPENDIX  D 

STABILISE  AS&L1SIS  FOE  THE  FINITE  DIFFERENCE 
ANALOG  OF  TBS  LINEAR  CONDUCTION  EQUATION 

A  difference  scheme,  such  as  that  obtained  by  assigning  values  to  r 
and  s  in  Etf,.  (25)  above,  will  be  termed  stable,  for  given  values  of 
the  space  and  time  increments  if,  initial,  truncation,  or  boundary  errors 
introduced  into  the  analysis  do  not  grow  without  bound  as  the  number  of 
time  steps  increases*  Thus  a  sufficient  condition  for  stability  iB  that 
the  ratio  of  the  total  error  at  the  end  of  any  computational  step  to 
that  at  the  beginning  does  not  exceed  one  in  modulus*  A  stability 
analysis  will  as  performed  for  Eq>  (25)  for  varying  r  and  s  to 
determine  which  values  provide  stable  schemes*  Stability  will  bo  further 
characterized  as  conditional  or  unconditional  according  to  whether  or  not 
it  depends  upon  the  boundedness  of  the  ratio  R  «= 

n 

$sing  B$.  (20/  above  and  its  radial  analog  for  d^=dB=dg=djj=h,  we  can 
rewrite  (25)  in  the  following  form  (the  indices  m  and  n  are 
to  indicate  spatial  subscripting  rather  than  i  and  j  to  avoid 

x confusion  below  with  the  imaginary  number  i  «=  /-l)t 

\, 

{^'4  *  -  {l^2sta-r)^(l-U)«J]}rmin,k  .  (B-l) 

If  we  igaoro  the  effect  of  boundary  conditions  and  lump  all  errors  intro¬ 
duced  into  the  analysis  at  any  point  into  a  3 ingle  variable, 


...•£■  ^  T';-'  ; 


8  ,  ,  then  we  see  that  *  „  .  will  satisfy  Bq.  (2-l)  within  terms 

o,n,k'  a»n,k 

of  order  h2*(rvs»l)AT  +  AT2*  Assuming  a  separated  solution  for  c 


m»u,k 


(see  Ref.  20,  for  example)  of  the  form 
~  ,  im»  h  in$  h 

e  =  V  A  e  **  e  ^  ,  (U-2) 

m,n,k  L  p,q  p,q  '  ' 

p.q 

where  p,q  take  on  integral  values,  then  if  we  can  show  that  (1^  ^|s  1 

for  all  p  and  q  this  will  imply  that  g  i,  This  is  the 

m,n,k 

sufficient  condition  stated  earlier  for  stability. 


Substituting  Eq.  (D-2)  into  (D-l),  equating  termB  of  like  indices  p 
and  q  and  eliminating  common  factors  we  obtain:  , 

{l-A[rt2«62]  }  ^Vnf\q=  {l^E[(l-r)4+(x-aj6=]}ei,,,^Vn^ 


Then  using  the  identities 


(D-J) 


h hi  eiBA  ->  -  *eiB0W  <f  , 

h2t2  olneh  =  .  te^sin2  •&  , 

El.  (B-3)  can  be  rewritten  in  the  following  form: 

l-4R[(l~r)sin2  &  +  (l-s)sin2 


1+4rJV  din2  ^  +  s  sin2 


(D-4) 


(D-5) 


U>-6) 


I 

mk  A 


We  wish  to  see  when  |7||  s  1  or  -1  sl]s  1.  Because  r,  s  a  1,  we 
see  that  T)  s  1  is  satisfied  for  all  permissible  values  of  r,s,  and 
Rj  T|  i  -1  is  satisfied  when 

R^(l-2r)sin2  +  (l-2s)sin2  £  £  a  (D-7) 

It  is  easy  to  see  that  (D-7)  holds  for  all  R  >  0  if 

r,  s  *•£•  •  (B-8) 


Thus,  if  r  and  s  satisfy  (D-8),  Eq.  (D-l)  determines  an  uncondition¬ 
ally  stable  soheme.  If  r,s  <  £,  then  (D-7)  holds  only  for 


Hs4iEfen 


(D-9) 


\  * 
and  we  haye  only  conditional  stability.  We  may  similarly  show  that  for 

one  of  the  parameters,  say  r,  less  than  %  and  a  £  -£•,  (D-7)  is 

satisfied  when  R  <  a01*  again  we  have  only  conditional 

stability.  These  cases  cover  all  possible  schemes  determined  by  Eq. 

(25)  for  which  r  and  s  are  invariant  and  we  therefore  see  that  the 

only  ones  which  satisfy  the  sufficient  condition  of  non-amplification 

for  all  values  of  R  are  those  given  by  (D-8)* 


There  is  another  class  of  schemes  in  which  values  for  r  and  s  are 
switched  after  every  time  step,  which  we  term  "alternating  schemes'?#  In 
particular,  we  will  consider  a  restricted  class  of  alternating  schemes 


HI 


for  which  s  =  1-r  ard  r  >  •£■  in  odd  time  steps.  The  alternating 
direction  method  tf  Peaceaan  and  Rachford  falls  within  this  class  and 
is  obtained  when  r  =  1  and  s  =  0  in  odd  time  steps  and  r  =  0  and 
s  =  1  in  even.  To  analyze  the  stability  of  this  restricted  class  of 
alternating  schemes,  we  investigate  the  following  two  expressions,  one 
for  an  odd  step  ard  one  for  an  even  which  are  analogous  to  Eq*  (D-6) t 


1-4rF (l-r  )sin2  ~  +  r  sin2 


^  =  f  *  (l-r0)sin2  &]  ' 


(]>-10) 


l-4s|ro8in2  +  (l-»ro)ain2 
^von  l+4R[(l-ro)sin2  +  rQsin2  &■'  * 


(D-ll) 


where  r  has  the  .value  rQ  in  odd  steps  and  l-rQ  in  even.  It  has 
already  been  demonstrated  that  neither  nor  !\venl  are 

bounded  by  1  for  all  values  of  R  >  0  but  only  for  values  of 
R  £  2X2^1)  s  2(l~2r)  *  r03P®ctively.  However,  it  can  readily  be 

shown  that  |\^\Ven*  5  *  *or  ®  >  provided  that  R  has  the 
same  value  in  both  Eq.  (D-IO)  and  Eq.  (D-ll),  i.3.f  provided  that  At 
be  held  constant  for  two  time  steps.  Therefore,  we  oan  assert  the 


unconditional  stability  of  all  such  restricted  alternating  schemes  for 
solving  Eq.  (D-l).  The  simplest  for  computational  purposes  and  there¬ 
fore  the  one  we  have  chosen  is  the  Peaceman-Rachford  method,  since  it 
is  the  only  one  which  generates  tridiagonal  system  of  equations. 


♦Since  the  writing  of  this  report  a  stability  analysis  was  performed 
covering  all  possible  alternating  schemes,  removing  the  restriction 
r  -t  s  =  1.  The  criterion  determined  for  unconditional  stability  of 
an  alternating  scheme  is  r  +  s  £  1  provided  that  each  value  of 
AT  employed  be  held  constant  for  at  least  two  time  steps.  It  can 
further  be  shown  that,  when  r  +  s  <  1,  the  condition  on  R  required 
for  stability  is  just  inequality  (D-9).  Clearly,  alternating  schemes 
provide  no  advantage  over  those  for  which  r  and  a  are  invariant, 
with  the  single  exception  of  the  alternating-direction  method  which  is 
computationally  more  convenient. 
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©Isl  repos  tuttahsr  by  which  the  <ocs*sat  will  bo  t&tltfied 
sod  co strolled  by  tha  erlglaatlag  activity.  This  tuaBhar  atuat 
be  unique  to  thla  report 

ca.  OTK33  REPORT  KUX3SR(3).  If  t!> .  <*jm  has  fesea 
nteignad  any  ether  rspert  number*  (either  bp  the  on'tntter 
at  kp  the  epeneet),  also  eater  thla  ourabarfa). 

19.  AVAIL ARU.irr/‘  ISTCATKJIt  JfeKCBO  Eater  any  It®. 
Uaiierta  oa  fyrth*'  divWmJaaUoa  of  the  report,  ether  then  those 


imposed  by  eternity  desslflcetloa,  using  standard  statantats 
such  as: 

(1)  "Qualified  requesters  may  Obtain  copiva  of  this 
report  bom  DDC." 

(2)  "PoraJgn  announcement  and  diatsa  action  at  thla 
report  by  D9C  Is  act  Authorised," 

(3)  "B.  &  Government  agendas  may  obtain  copies  of 
this  report  directly  from  CSC.  &het  qualified  OBC 
ustTs  shall  request  through 


*'U.  &  military  agaaclss  way  obtain  eopita  of  thla 
report  directly  from  HOC.  Othav  Qualified  ogvrg 
ahall  raquMt  through 


(3)  "All  dletrCcaUon  of  title  report  la  eoatfollod  Qud- 
If  ltd  ODC  uf  am  shall  roquoat  through 

_ _ _ _ _ _ 

If  tho  report  baa  boao  furnlahod  to  tho  Offlcv  of  Tochalcal 
* arvicos,  Oepartoeat  of  Conatcrca,  for  aalo  to  the  public,  ladh 
cate  this  fact  and  aster  tho  price,  if  known.  . 

It  tUP^EiEfTA^  rxnrta  Use  for  additional  avion* 
tory  cotes. 

U  C’OJtSCRItfQ  mUTARY  ACTIV2TT:  Enter  the  na»e  of 
the  dcpertasoatel  project  office  rr  laboretcry  eponaorto*  fpay- 
ieg  tot)  tho  reaaauch  sad  develops®*.  Isdstdo  address, 

13-  AC3T3ACT:  Eater  ea  abstract  giving  a  brief  sad  factual 
suauaary  «f  the  document  Indicative  of  tho  report,  aver,  though 
It  may  alto  eppoar  alsowfeora  la  the  body  of  the  techalcalro* 
port  tf  eddltleae!  opace  la  legalni  a  coatlnuatloa  aheet  ahall 
bo  attached. 

It  la  hifhly  dealrahle  that  tho  ahatract  of  claeelflod  reports 
be  unclaaslfled.  Each  paragrepb  of  the  ahatract  ahall  end  with 
aa  Indication  of  fha  raliltary  security  classification  of  tbs  in¬ 
formation  la  the  paK’graph,  rsprssaatsd  aa  tTSj.  1 3j.  (Cj.  «r  ruj 

There  la  so  Eatltettoa  on  the  length  of  the  abstract,  flow- 
over,  the  suggested  leo$5«  la  from  12.3  to  229  worda. 

14  ESY  WKDSr  Sc?  wok1«  are  t*ihiUc»Uy  meamnsfal  tens* 
abort  phnt8t*<  that  el waeP*n»*  «  report  and  mey  he  used  oa 
lade  aatrles  for  cats  logins  the  report.  Key  worda  tasat  bo 
oelectel  so  that  no  eecurlty  ciaeslflcauoo  ta  required.  Ideatl-. 
Car.'  eu  th  ea  sqalpanoBt  taodsl  d#tI»*tlon,  teeCe  aeme,  military 
project  redo  rtaate,  neocraphic  loeetloa,  oar/  b *  art  ea  hsy 
words  hut  will  ho  fcuewsd  by  sa  ladlcatlen  of  techsicsi  con- 
tec*  *3tf  eeslgscwat  ot  Uohe.  rales,  and  mights  is  opunoal.  i 


Socsrity  ClsselOccuca 


